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ABSTRACT 

A  fully  factorized  two-dimensional  Navier-Stokes  flow  solver  has  been 
developed  and  applied  to  the  problem  of  predicting  subsonic  airfoil  flutter  in  the 
light  stall  regime.  The  inviscid  fluxes  are  evaluated  with  a  central  difference 
ADI  scheme  and  fourth  and  second  order  numerical  dissipation  is  used  to  obtain 
oscillation-free  solutions.  The  performance  of  algebraic  and  one-equation 
turbulence  models  in  predicting  separated  flow  is  explored  for  computing  high 
Reynolds  number  steady  flow  and  unsteady  flows  over  an  oscillating  NACA  0012 
airfoil.  Comparisons  of  the  computed  results  with  available  experimental  data 
indicate  that  even  though  the  lift  response  is  fairly  well  predicted,  the 
computation  of  the  pitching  moment  hysteresis  loops  is  very  sensitive  to 
tuibulence  modeling.  Results  computed  with  several  current  models  are  in  good 
agreement  whenever  the  steady  stall  angle  is  exceeded  only  slightly.  However, 
they  fail  to  capture  the  vortex  shedding  process  leading  to  the  onset  of  stall 
flutter. 
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I.    INTRODUCTION 

A.  BACKGROUND 

The  development  of  numerical  solution  methods  for  two-dimensional 
Navier-Stokes  equations  during  the  past  few  years  provides  new  tools  for  the 
investigation  and  prediction  of  airfoil  flows.  Of  great  interest  is  the  study  of  flow 
separation  on  airfoils  in  unsteady  motion  which  is  usually  referred  to  as  the 
dynamic  stall  phenomenon.  McCroskey  et  al  [Ref.  1]  performed  a  series  of 
careful  experiments  which  serve  as  valuable  benchmark  data.  More  recently, 
Lorber  and  Carta  [Ref.  2]  contributed  important  additional  experimental 
dynamic  stall  information  for  a  Sikorsky  airfoil.  They  also  investigated  incipient 
torsional  stall  flutter  [Ref.  3]  and  found  that  small-amplitude  airfoil  oscillations 
near  static  stall  may  be  unstable. 

For  a  pitching  airfoil  the  instantaneous  work  done  on  the  fluid  by  the  airfoil 
due  to  its  motion  is  the  product  of  the  pitching  moment  about  the  axis  of  rotation 
and  the  differential  change  in  angle  of  attack.  This  product  usually  is  a  positive 
quantity.  However,  if  the  net  work  per  cycle  of  oscillation  were  to  become 
negative  then  the  fluid  would  be  doing  work  on  the  airfoil.  Once  the  airfoil 
begins  to  extract  energy  from  the  freestream,  the  amplitude  of  the  oscillation 
will  grow  and  finally  diverge.  This  condition  is  known  as  stall  flutter.  Normally, 
flutter  of  an  airfoil  is  due  to  a  combination  of  torsion  and  bending.  However,  in 
this  case  flutter  is  caused  by  a  single  degree  of  freedom  oscillatory  motion. 
During  a  cycle  of  oscillation  the  lift  coefficient  and  the  pitching  moment 
coefficient  plotted  versus  angle  of  attack  produce  a  hysteresis  loop.    It  is  the 


pitching  moment  hysteresis  loop  that  provides  an  indication  of  incipient  stall 
flutter.  As  shown  by  Carta  and  Niebanck  [Ref.  4],  clockwise  pitching  moment 
loops  represent  negative  aerodynamic  damping  and  therefore  cause  oscillations 
of  a  free  airfoil  to  grow  in  amplitude,  while  counterclockwise  loops  cause  such 
oscillations  to  decay.  Hence  torsional  flutter  will  occur  as  soon  as  the  area  of 
the  clockwise  loop  exceeds  that  of  the  counterclockwise  loop.  The 
aerodynamics  of  small-amplitude  airfoil  torsional  oscillations  near  stall  need  to 
be  investigated  experimentally  and  computationally  in  order  to  examine  the  stall 
flutter  mechanism. 

B.    PURPOSE 

The  first  objective  of  this  investigation  is  to  test  an  unsteady,  compressible 
Navier-Stokes  code  (Ns2.f)  using  an  Alternating-Direction-Implicit  (ADIj 
scheme  based  on  the  the  Beam-Warming  [Ref.  5]  approximate  factorization 
method  to  determine  its  ability  to  obtain  realistic  airfoil  flow  solutions  for  a 
variety  of  flow  regimes  is  examined.  The  accuracy  of  the  numerical  solution  is 
investigated  by  comparing  the  computed  solutions  with  available  experimental 
data.  Test  cases  include  steady-state  flow  solutions  at  various  flow  speeds  and 
angles  of  attack,  as  well  as,  unsteady  flow  solutions  over  rapidly  pitching  and 
harmonically  oscillating  airfoils.  In  addition,  the  accuracy  of  the  Navier-Stokes 
solutions  are  further  explored  by  comparing  each  case  with  an  unsteady, 
inviscid,  incompressible,  panel  code  (U2diif.f),  and  with  a  steady, 
incompressible,  viscous/inviscid  interaction  code  (Incompbl.f).  The  final  and 
principal  objective  of  this  study  is  to  determine  the  influence  of  mildly  separated 
flow  during  part  of  a  harmonic  oscillation  cycle  on  blade  stability.  Results  are 
presented  for  NACA  0012  airfoils  showing  the  influence  of  various  parameters. 


Row  field  details  are  also  included  in  order  to  provide  a  better  understanding  of 
certain  flow  features  that  may  lead  to  stall  flutter. 


II.  GOVERNING    EQUATIONS 

The  continuity  equation,  the  momentum  equation,  and  the  energy  equation 
must  be  solved  simultaneously  in  order  to  obtain  a  flow  solution  of  a 
compressible  viscous  fluid  about  a  body.  A  complete  derivation  of  these 
equation  can  be  found  in  various  texts,  eg.,  Anderson  [Ref.  6].  However,  main 
steps  of  the  derivations  are  presented  in  the  following  sections. 

A.   CONTINUITY    EQUATION 

The  continuity  equation  is  the  result  of  applying  the  physical  principle  of  the 
conservation  of  mass  to  a  finite  control  volume  fixed  in  space.  Simply  stated 
the  net  flow  out  of  a  control  volume  through  its  bounding  surface  must  be  equal 
to  the  time  rate  of  change  of  the  mass  inside  the  control  volume.  For  any 
arbitrary  control  volume,  the  continuity  equation  can  be  expressed  as 


§?+v(Pv)=o 

which  for  a  two-dimension  Cartesian  Coordinate  system  becomes 
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B.    MOMENTUM    EQUATIONS 

The  equation  for  the  conservation  of  momentum  is  obtained  by  applying 
Newton's  second  law,  which  state  that  the  net  force  acting  on  a  fluid  particle  is 
equal  to  the  time  rate  of  change  of  linear  momentum  of  the  fluid  particle.    In 


Cartesian  coordinates  the  momentum  equations  can  be  expressed  as  follows 
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and  in  the  y-direction. 
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Stress  terms  are  as  follows 
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C.    ENERGY    EQUATION 

The  conservation  law  form  of  the  energy  equation  is  derived  by  applying 
the  first  law  of  thermodynamics  (dE=dQ+d\V)  to  a  fluid  particle.   This  leads  to 


§  +  |[(E+p)u]  +  |[(E  +  p)w]  = 
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where  the  total  energy  per  unit  volume  and  the  heat  flux  are  given  by 


(6) 
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Here  k  is  the  thermal  conductivity. 


D.    CONSERVATION  FORM   OF  GOVERNING   EQUATIONS 

The  starting  point  for  the  numerical  algorithm  which  is  presented  in  the  next 
section  is  the  strong  conservation  law  form  of  the  two-dimensional  Navier- 
Stokes  equations.  The  non-dimensionalized  vector  form  of  the  governing 
equations  in  conservation  law  form  for  a  Cartesian  coordinate  system  is: 


where, 
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Here  the  pressure  is  related  to  the  conservation  variables  of  q  by  the  equation 


p  =  (7-l)[e-0.5p(u2  +  v2)] 


(11) 


where  y=  1.4  is  the  ratio  of  specific  heats,  and  a  is  the  local  speed  of  sound 
given  by,  a^  =  yp/  p. 

The  density  and  the  velocities  are  non-dimensionalized  with  the  freestream 
density  poo  and  the  free  stream  speed  of  sound  aoo,  respectively.  The  total 
energy  is  normalized  by  a^pcc-   The  time  (t)  scales  as  t*  =  t  aoo  /  c,  where  (c) 

is  a  characteristic  length,  such  as  the  chord  length.  The  Euler  equations  are 
obtained  from  Equation  (8)  by  dropping  the  viscous  terms.  The  strong  form  is 
chosen  because  it  enables  shock  capturing. 


III.     SOLUTION   METHODS 

A.   POTENTIAL  FLOW  METHOD  (U2DIIF.F) 

For  inviscid  incompressible  flow  the  conservation  equations  reduce  to  the 
Laplace  equation  which  can  be  solved  in  a  number  of  ways.  We  present  here  a 
brief  outline  of  the  widely  used  panel  method  for  both  steady  and  unsteady 
airfoil  flow. 

1.      Panel  Method  for  Steady  Airfoil  Flows 

The  frame  of  reference  for  the  formulation  of  the  steady  flow  problem 
is  a  fixed  (x,y)  coordinate  system  located  at  the  leading  edge  of  the  airfoil.  The 
freestream  velocity  is  represented  by  +V©o  and  the  x-axis  passes  from  the 
leading  edge  through  the  trailing  edge  of  the  airfoil. 

In  order  to  formulate  a  method  for  computing  the  flow  around  an 
airfoil,  the  airfoil  surface  is  divided  into  (n)  straight  line  segments  or  panels.  The 
end  points  of  the  panels  are  called  nodes,  and  since  there  are  (n)  panels  there 
are  then  (n+1)  nodes.  Complex  airfoil  geometries  can  therefore  be  modeled 
with  a  greater  number  of  nodes  and  panels.  The  trailing  edge  panel  is  the  first 
panel,  and  the  next  panel  continues  on  the  lower  surface  in  a  clockwise  manner 
until  the  nth  panel  is  reached  at  the  trailing  edge  of  the  upper  surface. 

Now  that  the  frame  of  reference  has  been  specified,  two  additional 
vectors   must   be   defined,   the   unit   normal   vector   nj,    which   is   always 
perpendicular  to  the  (im  )  panel  and  directed  outward  from  the  airfoil  surface 
and  the  unit  tangential  vector  ti,    which  is  parallel  to  the  (itn)  panel  and  is 
directed  from  the  (n)  node  to  the  (n+1)  node. 


Two  types  of  singularities  are  sufficient  to  model  lifting  airfoil  flows. 
U2diif.f  places  uniform  source  distributions  qj  and  a  uniform  vorticity  distribution 
y  on  each  of  the  j-panels.  Both  singularities  satisfy  Laplace's  equation,  a  linear 
homogeneous  second  order  partial  differential  equation.  Since  the  solutions  to 
linear  PDE's  can  be  superimposed  on  each  other,  a  simple  flow  can  be  added  to 
another  simple  flow  and  so  on,  until  a  very  complicated  flow  field  is  created. 

The  flow  field  around  an  airfoil,  represented  by  the  velocity  potential, 
can  be  constructed  from  the  potential  of  the  freestream  flow  added  to  the 
velocity  potential  of  source  and  vorticity  distributions,  hence 
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where, 


%  =  V00+(x  cosa  +  y  sin  a) 
summation  of  the  three  potentials  yields 
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Equation  (16)  can  now  be  evaluated  at  any  point  in  the  flow  field  by  evaluating 
the  integrals  along  the  airfoil  contour  (s),  where  the  flow  field  point  is  located  at 


a  distance  (r)  at  an  angle  (0)  measured  from  a  point  on  the  airfoil.    Introducing 
n  panels  and  then  summing  the  contributions  of  each  panel  yields 
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Once  the  boundary  conditions  are  defined  the  system  of  (n)  equations  can  be 
solved  for  (n+1)  unknowns. 

Next,  it  is  useful  to  introduce  the  concept  of  influence  coefficients.  An 
influence  coefficient  is  the  velocity  induced  at  a  point  in  the  flow  field  (field 
point)  by  a  unit  strength  singularity,  source  or  vorticy,  placed  anywhere  within 
this  field.  U2diif.f  places  these  points  on  each  panel.  A  detailed  description  of 
the  use  of  influence  coefficients  is  found  in  Teng  [Ref.  7]. 

For  the  steady  flow  problem  the  influence  coefficients  are  given  in 
Table  1. 

TABLE    1.   INFLUENCE    COEFFICIENTS 


Any 

Normal  component  induced  at  i1*1  control  point  by 
unit  source  distribution  on  jtn  panel. 

Atjj 

Tangential  component  induced  at  itn  control  point 
by  unit  source  distribution  on  jth  panel. 

B"ij 

Normal  component  induced  at  itJl  control  point  by 
unit  vorticity  distribution  on  jth  panel. 

Btjj 

Tangential  component  induced  at  itn  control  point 
by  unit  vorticity  distribution  on  jtn  panel. 

6ij 

Angle  of  control  point  (i)  and  panel  (j)  measured 
from  x-axis. 
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There  are  only  two  boundary  conditions  that  must  be  satisfied  in  order 
to  solve  the  lifting  airfoil  flow  problem.  The  first  boundary  condition  is  the 
condition  that  the  flow  remains  tangent  to  the  airfoil  surface  and  the  second 
condition  is  that  the  pressures  on  the  upper  and  lower  trailing  edge  panels  must 
be  equal.  This  condition  is  known  as  the  Kutta  condition.  Pressure  is  related  to 
velocity  through  Bernoulli's  equation  for  steady  potential  flow,  which  allows  the 
Kutta  condition  to  be  expressed  as 

/y  tangent^  __[ytangent) 

^lower  'upper  M8^ 

The  flow  tangency  condition  is  expressed  as 

(ynormal)    =  q         i  =  ^2  n 

x  'i 

(19) 
Equation  (18),  the  Kutta  condition,  using  the  influence  coefficient  concept  is 
expressed  as  follows 

-X[A1tjqJ-rZ[B;j]-V„cos(a-01)  = 

j=i  j-i 

XkflJ+rXkJ  +  Vocos(a-0j 

j=i  H  (20) 

Equation  (19)  the  flow  tangency  condition,  becomes,  using  the  influence 

coefficient  form 


S[A!,jqj]+rZ[B;J+V_sin(a-01)  =  a     i=l,2,...,n 

H  H  (21) 
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These  equations  can  then  be  written  in   the  following  matrix  form 
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and  solved  with  the  method  of  Gauss  Elimination  with  Partial  Pivoting. 
2.      Unsteady   Numerical   Formulation 

A  brief  and  concise  description  of  the  unsteady  numerical  formulation 
is  found  in  Krainer  [Ref.  8].  To  model  unsteady  flow  around  an  airfoil  using  the 
panel  method,  N  unknown  source  strengths  and  one  unknown  vorticity  strength 
are  located  along  the  airfoil,  a  wake  panel  of  unknown  vorticity  strength,  length 
and  orientation  is  attached  to  the  trailing  edge  of  the  airfoil.  This  makes  a  total 
of  N+4  unknowns.  The  flow  tangency  requirement  at  each  of  the  N  panels 
represents  a  system  of  N  equations.  Using  influence  coefficients  and  using  the 
subscript  k  to  count  time,  the  tangential  flow  condition  of  the  i1*1  element  is 
written  as  follows 


I  AiM  +rk£Br.+  (V.  +  fUWi  +  VWjJ  +  QdKxi-xiftn. 

j=lL  KJ       j=l        L  l 

+(7w)kte+1)  +  S[te)(r_1  - rjj  =  0    i  =  1,2,..., n 


m=l 


(23) 


where  Voo  is  the  mean  velocity,  (U(t)i  +  V(t)j)  is  a  time  dependent  translational 
velocity  and  Q(t)  is  a  rotational  velocity.  The  vectors  i  and  j  are  in  the  airfoil 
fixed  coordinate  system. 
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The  Kutta  condition  is  a  single  equation  that  determines  the  circulation 
around  the  airfoil.  In  this  case  the  pressures  at  the  upper  and  lower  surface  are 
equated  at  the  midpoints  of  panels  adjacent  to  the  trailing  edge  (i.e.,  the  first  and 
last  panel) 


[(vl)kf  +  [(vi)kf=2 


2_J"3(ft,-ft)l  _„r<?rl  _„,yk-rk-, 

(24) 
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dt 
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The  Helmholtz  theorem  is  another  single  equation  that  provides  a 
relation  for  the  strength  of  the  vorticity  distribution  along  the  wake  element. 
The  Helmholtz  theorem  states  that  any  change  in  circulation  around  an  airfoil 
must  be  countered  by  a  change  in  vorticity  in  the  wake  of  equal  magnitude  but 
of  opposite  sign.   This  theorem  can  be  written  as  follows 

Ak(7w)k  =  /(yk-i-7k) 

(25) 

where  /  is  the  length  of  each  individual  wake  element  and  y  its  vorticity 
strength.  Therefore,  the  vorticity  in  the  wake  element  is  equal  to  the  negative 
change  in  circulation  around  the  airfoil  with  respect  to  the  k-1  timestep. 

To  acquire  the  two  additional  equations  required  to  solve  for  the  N+4 
unknowns,  assumptions  about  the  geometry  of  the  wake  panel  are  made.  First, 
the  wake  panel  is  oriented  in  the  direction  of  the  local  resultant  velocity  at  its 
midpoint,  as  viewed  in  the  (moving)  airfoil  frame  of  reference. 


(vy) 


'^k  (26) 
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(vyw)k  and  (vxw)k  are  the  x-  and  y-velocity  components  at  the  midpoint  of 
each  wake  panel.  Second,  the  length  of  the  wake  panel  is  assumed 
proportional  to  the  magnitude  of  the  local  resultant  velocity  at  its  midpoint  and 
to  the  timestep. 

Ak=^(v:)J2+[(v:)k]2(tk-tk.1) 

(27) 
The  nonlinearities  in  the  Kutta  condition  and  in  the  wake  panel 
assumptions  necessitate  an  iterative  solution  procedure. 

B.    VISCOUS/INVISCID    INTERACTION    METHOD 
(INCOMPBL.F) 

The  theory  and  numerical  methods  presented  in  this  section  are  taken  from 
the  work  of  Cebeci  and  Bradshaw  [Ref.  9]  and  the  investigations  performed  by 
Krainer  [Ref.  8]  and  Snir  [Ref.  10]  As  stated  in  the  introduction  the  results 
from  the  Navier-Stokes  code  (Ns2.f)  were  compared  with  a  viscous/inviscid 
interaction  method  code  (Incompbl.f)  made  available  by  Cebeci.  An 
abbreviated  description  of  this  method  is  given  in  this  section.  For  a  complete 
description  see  the  original  publications,  References  8  and  9. 

The  governing  principle  behind  the  viscous/inviscid  interaction  method  is  an 
approximation  to  the  Navier-Stokes  equation  that  allows  a  flowfield  to  be 
divided  into  an  inner  viscous  region  and  an  outer  inviscid  region  when  the 
Reynolds  number  is  sufficiently  large.  This  concept  of  a  thin  viscous  layer  near 
the  surface  of  a  body  and  an  outer  region  where  these  viscous  effects  are  small 
compared  to  the  inertial  effects  is  known  as  the  boundary  layer  theory.  Unlike 
the  Laplace  equation,  which  is  the  governing  equation  for  potential  flow,  the 
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thin  shear  layer  equations  are  nonlinear 


eta     c^v 
dx     dy 


eki        dv  dup         d  I7t     v.^ctal 


(28) 
with  the  boundary  conditions 

y=0,  u=v=0, 

y=oo,  u=ue(x)  (29) 

The  direct  boundary  layer  method  solution  to  the  boundary  layer  equations 
by  the  Keller  Box  Method  involves  four  steps.  First  the  boundary  layer 
equations  are  transformed  into  a  system  of  first  order  differential  equations. 
Next  the  Keller  box  method  is  used  to  approximate  the  first  order  differential 
equations  by  simple  centered  differences  and  two-point  averages,  using  values 
at  the  corners  of  one  difference  molecule  only.  Newton's  method  is  used  to 
linearize  the  resulting  algebraic  equation.  The  Keller  block  elimination  method 
is  then  used  to  solve  the  resulting  block  tridiagonal  system.  The  procedure 
described  above  does  not  provide  solutions  to  flows  that  are  separated  or  have 
regions  of  reverse  flow. 

1.      Interaction   Method 

The  interactive  boundary  layer  method  provides  a  special  coupling 
between  the  inner  viscous  flow  and  the  outer  inviscid  flow,  which  enables 
reverse  and  separated  flows  to  be  calculated.     In  such  areas,  the  external 
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velocity  is  substantially  changed  by  the  viscous  effects  and  can  no  longer  be 
considered  as  a  known  boundary  condition  for  the  boundary  layer  flow. 

The  general  approach  to  the  solution  is  the  same  as  for  the  direct 
method  with  modifications.  Since  the  outer  flow  is  unknown,  the  velocity  at  the 
edge  of  the  boundary  layer  is  given  by  the  interaction  law  and  is  written  as 

71     d<>  x_(=  (30) 

where  ue(x,ye)  is  the  total  velocity  at  the  edge  of  the  boundary  layer,  uei(x)  is 
the  velocity  computed  by  the  inviscid  panel  method,  5*  is  the  displacement 
thickness,  and  the  integral  term  is  known  as  the  Hilbert  integral. 

The  following  transformations  are  used  to  transform  the  boundary 
layer  equation  into  a  system  of  first  order  differential  equations 


x  I —  y  u0L 


c/      \         1        ,/       x  ue(x,y) 

f(x,77)  =  -==!//(  x,y),       w  = 

Vu0vx  u0 

The  boundary  layer  equation  takes  the  form 

f=U,       U'=V,       W'=0, 
with  boundary  conditions: 


(31) 


(32) 


n=o,  u(x,o)=o,  f(x,o)=o, 

Tl=Tie,  U(x,Tle)=W(x,Tle)  (33) 
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The  velocity  at  the  edge  of  the  boundary  layer  now  becomes 

u0       n     d^\u0  )x-$  (34) 

The  finite  difference  box  method  is  used  to  solve  the  equations,  in  the 
same  way  it  was  used  for  the  direct  case,  but  with  two  additions.  First  in  areas 
of  flow  reversal  the  u3u/3x  is  omitted  to  assure  stable  integration.  And  second, 
the  edge  velocity  is  approximated  by  the  relation  presented  in  the  next  section. 

By  using  central  differencing  to  approximate  the  differential  equations, 
a  system  of  nonlinear  algebraic  equations  is  obtained  for  the  unknown  variables 
(which  are  f,  U,  V,  and  W).  To  solve  the  system  of  equations,  the  system  is 
linearized  by  the  Newton  iterative  procedure,  and  the  resulting  linear  system  is 
solved  for  the  new  unknown  variables  which  are  the  increments  5f,  SU,  5V,  and 
6W. 

The  solution  procedure  is  repeated  until  the  change  in  the  increments 
is  negligible  compared  to  the  preceding  iteration.     The  iterative  process  is 
performed  again  at  the  next  downstream  station. 
2.     Interactive   Model 

The  interactive  model  is  used  to  couple  the  boundary  layer  to  the 
external  flow.  It  is  needed  in  areas  where  strong  interaction  occurs,  and  both 
the  boundary  layer  and  the  outer  flow  must  be  solved  simultaneously.  The 
external  velocity  is  assumed  to  consist  of  a  potential  flow  term  and  a  correction 
term  due  to  viscous  effects.  [ue(x)=uel(x)+ue5(x)] 
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The  normal  velocities  at  the  surface  of  the  airfoil,  induced  by  sources 
distributed  on  the  surface,  displace  the  streamlines  from  the  surface  in  the  same 
way  that  the  actual  boundary  layer  displaces  them. 

Several  assumptions  are  made  in  order  to  express  the  correction  term 
in  the  form  of  the  Hilbert  integral.  First,  the  surface  is  approximated  to  be  a  flat 
plate,  and  the  normal  velocity  is  then  half  the  local  source  strength  o(x). 
Second,  the  inviscid  velocity  does  not  change  across  the  boundary  layer. 
Therefore,  the  local  horizontal  velocity  induced  by  the  source  distribution,  is  the 
correction  term  to  the  inviscid  velocity,  and  can  be  represented  by  the  Hilbert 
integral 

Ij£(|^aj*(u/)-iL 

/rJ  x-£  *     k>  d£v        'x-§  (35) 

The  integration  is  carried  out  on  all  the  sources  on  the  surface.    The  Hilbert 
integral  is  then  approximated  by  a  finite  series 


y^^-fr^r 


(36) 

where  cjk  is  a  matrix  of  interaction  coefficients  which  are  functions  of  the 
geometry  only. 

Since  the  computation  of  ue5  involves  values  of  5*  downstream  of  the 
current  x  location,  which  are  not  known  yet,  these  terms  are  taken  from  the 
previous  iteration  using  a  relaxation  formula. 
3.     Turbulence   Model 

The  turbulence  model  used  in  the  code  is  the  Cebeci-Smith  model 
[Ref.  11]  and  Michel's  method  is  used  to  predict  the  transition  from  laminar  to 
turbulent  flow. 
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C.    EULER   AND   NAVIER-STOKES   CODE   (NS2.F) 

In  order  to  facilitate  the  numerical  implementation  for  arbitrary  complex 
flow  domains,  the  Navier-Stokes  equations  are  transformed  to  a  generalized 
coordinate  system,  (£,Q,  using  the  transformations 

£  =  «k,z) 

C  =  C(x,z)  (37) 

The  grid  spacing  in  the  curvilinear  space  is  uniform  and  of  unit  length  so  that 
unweighted  differencing  schemes  can  be  employed  for  the  numerical 
implementation.  The  grid  points  in  the  Cartesian  system,  referred  to  as  the 
physical  domain,  correspond  via  a  one  to  one  relationship  to  the  points  of  the 
curvilinear  transformed  system,  referred  to  as  the  computational  domain. 
Singularities  of  the  transformation  may  occur  on  the  computational  boundaries. 
Such  singularities  occur  for  grids  with  multiple  connected  regions.  The 
transformation  of  the  governing  equations  from  the  physical  domain  to  the 
computational  domain  is  obtained  in  most  cases  by  numerically  evaluating  the 
metrics  and  the  Jacobian  of  the  transformation.  The  derivatives  with  respect  to 
the  x,  y  variables  can  be  expressed  in  terms  of  the  new  variables  by  the  chain 
rule. 


a       j.    a      y,    d 

a       „    d      r    d 

*      "^       ^  (38) 
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The  metrics  and  the  Jacobian  of  the  transformation  are  given  by 


£x  =  J  Z£, 

5z  =  -J  x£, 


Cx  =  -J  z£, 
Cx  =  J  xg, 


J  = 


x^zc     x?z^ 


(39) 
(40) 

(41) 


The  Navier-Stokes  equations  written  in  generalized  curvilinear  coordinates 
retain  the  strong  conservation  law  form  expressed  by  Equation  (8).  The 
governing  equation  for  generalized  curvilinear  coordinates  is 


^L    it    £*i     JJ^El    £*L 
d  t +  %  +  <9f  "  Re    a*  +  d£ 


(42) 


The  conservation  variable  vector  Q  and  the  inviscid  fluxes  F  and  G  in  the 
transformed  space  are  given  by 


\P    1 


Tpu 


q"jlpwl'         =JpwU  +  ^p 


1  Tpw 

^lp\Vu+Cxp 

!'  g"j!pww+czP 


1 


L(e+p)U-£pj  |_(e  +  p)W-CtPj  (43) 

where  U  and  W  are  the  contravariant  velocity  components  along  the  £,  and  £ 
directions  respectively,  given  by: 


The  viscous  flux  terms  are  transformed  as 


Fv  =  j&Fv  +  £Gv),       Gv=-J-(CxFv  +  fzGv) 


(44) 


(45) 
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with  the  stress  terms  from  Equation  (10)  expressed  in  terms  of  the  transformed 
variables  £,  £  as 


^  =  3^(^+£uJ--(^+CzwJj 
^z=m[(£^+CzuJ+(^+CxwJ] 


£=UT„  +  WT     + 


^ 


g4=urxz  +  wTzz  + 


z    Pr(y-l) 


;.   ^a  eta 

„  da.       „  da. 

^-TF+Cz 


(46) 


Pr(7-1)L"Z  <?£     bz  (9C. 

When  the  thin  layer  approximation  [Ref.  12]  to  the  two-dimensional 
conservation  law  form  of  the  governing  equations  is  applied,  they  take  the 
following  form 


where, 


£§.  ^£  ^2  jJ^L" 

<9t+^  +  <9£  ~ReL<?C  J 


(47) 


rp  i 

«     llpu  I 
q=7lpwl' 


rpu         i 

jlpuU  +  &p       I 

jjpwu+4p     ' 

[(e  +  p)U-£pJ 


TpW  1 

JpWu+C,P 
G-jjpwW  +  CzP      ] 
|_(e  +  p)W-£pj 


(48) 
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The  viscous  flux  term  is  transformed  as 


{         °  ) 

g       J_    ^mlU?  +  (^/3)m2^ 
Jv     J    fi  m,wc  +  (/i/3)m2  fx 
^  m1m3  +  (/2/3)m2m4> 


(49) 


Here, 

1  <?  /  2      2 Wo  1       ^l 

ms  =  _^(n+w)/2+_(_j^j 

m4  =  C,u  +  ?2w  (5Q) 

1.      Numerical   Grids 

In  order  to  compute  flow  solutions  of  partial  differential  equations 
(such  as  the  Navier-Stokes  Equations)  with  finite  differences,  a  discretized 
version  of  the  physical  domain  must  be  generated.  During  the  course  of  this 
investigation,  several  different  NACA  0012  airfoil  C-type  grids  were  generated 
using  two  different  grid  generation  methods.  A  grid  generation  code  call 
GRAPE  that  uses  the  Poisson  differential  equation  was  used  in  the  early  stages 
of  investigation.  The  inviscid  and  viscous  grids  shown  in  Figures  1  and  2  were 
generated  from  an  algebraic  grid  generation  code.  Note  the  finer  resolution  of 
the  viscous  grid  near  the  airfoil  surface.  This  resolution  is  required  to 
accurately  represent  the  boundary  layer.  The  flow  solutions  presented  in 
Chapter  IV  were  computed  by  using  the  grids  shown  in  Figures  1  and  2.    The 
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advantage  to  using  an  algebraic  grid  ,  vice  a  PDE  method  or  a  conformal 
mapping  method,  is  that  it  is  computationally  efficient. 

A  simple  description  of  the  algebraic  grid  generation  method  is  offered. 
First  the  airfoil  surface  coordinates  are  unwrapped  to  form  a  simple  curve,  by 
separating  the  lower  trailing  edge  from  the  upper  trailing  edge.  This  curve  is 
the  starting  point  for  the  generation  of  the  computational  domain.  Lines  are 
drawn  normal  to  the  curve  and  grid  points  are  then  located  at  intervals  on  these 
lines.  Lines  may  be  clustered  near  the  leading  or  trailing  edge  of  the  airfoil;  or 
at  some  other  area  of  interest,  eg.,  a  shock  location.  The  resolution  of  the 
points  normal  to  the  surface  is  achieved  by  the  use  of  an  algebraic  function  that 
provides  uniform  stretching  normal  to  the  body  surface.  Once  the  grid  is 
created  in  the  computation  domain  it  is  then  transformed  back  into  the  physical 
domain  using  inverse  transformations. 
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Figure  1.     NACA  0012  Inviscid  151x64  C-Grid 
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Figure  2.     NACA  0012  Viscous  161x64  C-Grid 
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2.      The   Numerical    Scheme 

The  integration  in  time  is  performed  implicitly  with  a  second  order 
accurate  scheme  of  the  form 


n_„n  +  l     _n     At  d  ( A  n\  ,  At  dqa  ,  J  At2  ^      3^ 


^'^""-q-jKj^t^ 


(51) 

Using  this  trapezoidal  rule  differencing  scheme  to  approximate  the  time 
derivative  of  the  conservative  dependent  variable  q  vector,  the  unknown 
conservative  variables  at  the  n+l  timestep  are  given  by 


n  +  l  ,  AtJfdqY     (dq 


Aq "=q"  Ti+y< 


The  spatial  derivatives  of  the  governing  equations  are  discretized  using  central 
differences,  and  the  right  hand  side  of  Equation  (52)  becomes 


*in  =~w+sfij +(¥+sfir] 


2H-«        1    >      ™        1    >     J  (53) 

The  nonlinear  terms  F  and  G  at  the  n+l  timestep  are  linearized  as  follows 

rip" 

Fn+1  =  Fn  +^-Aqn+1  +0(At2)  =  Fn  +  AnAqn+1  +0(At2) 

<*!  (54) 

where  An  is  the  flux  Jacobian  matrix  given  in  Appendix  A.    Substitution  of  the 

linearized  form  of  the  flux  vectors  in  Equation  (54)  yields 


L  J  (55) 
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The  two-dimensional  operator  on  the  left  hand  side  of  Equation   (55)  is 
approximately  factorized  as  follows: 

(56) 
On  the  right  hand  side  of  Equation  (56)  an  explicit  dissipation  term  Dexpl,  is 
added  for  numerical  stability.  In  order  to  obtain  an  oscillation  free  solution  a 
second  order,  dissipation  term  Dimpi  is  added  to  the  left  hand  side  spatial 
operators.  Thus,  the  complete  discretized  form  of  the  governing  equations 
becomes 


This  equation  is  solved  by  performing  two  sweeps  as  follows: 


(57) 


(58) 


(59) 


{uj+f^B-k+(DiBP1)J}=^.r 


(60) 

During  each  sweep  the  following  linear  block  tridiagonal  system  of  equations  is 
solved 

^-i.k^-i.k+b^Aq.k+c^^kAq"^^^  (61) 
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where  Aq  =  Aq*n+1  and  f  =  (RHS)n  when  the  £,  direction  sweep  is  performed, 
and  Aq  =  Aqn+1  and  f  =  (RHS)*n  when  the  C,  direction  sweep  is  performed. 
The  block  matrices  have  the  form 


At  a        .a.  J& 

a1,k=yAlk-eimplAt— 


i+l,k 


b1>k  =  (l-h2f:impIAt)[l] 

At  -  Jlk 

Ci,k  =~T^i,k  "^mpA1"}- 
z  Ji-l,i 


(62) 
(63) 
(64) 


3.      Boundary  Conditions 

All  flows  were  computed  at  subsonic  free  stream  speeds.  For 
subsonic  inflow  and  outflow  boundaries  the  flow  variables  are  evaluated  using 
zero  order  Riemann  invariant  extrapolation.  At  the  inflow  boundary  there  is  one 
incoming  and  three  outgoing  characteristics,  therefore,  three  variables,  density 
(p),  normal  velocity  (wi),  and  pressure  (p)  are  specified  and  the  fourth  variable 
axial  velocity  (ui),  is  extrapolated  from  the  interior.  The  inflow  boundary 
conditions  are  given  by 


\   i»< 


Ul  =  (Rf  +  R2-)/2 


s,  =(p„/pj, 


ai=ilzH(R;-R-) 


wt  =w0 


(65) 
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where  R+i,  R~2  are  the  incoming  and  outgoing  Riemann  invariants  given  by 

Rj=.u.  +  2a-/(y-l)i       R"  =u2-2a2/(y-l) 

(66) 

At  the  outflow  boundary  there  is  one  incoming  and  three  outgoing 

characteristics  and  only  one  quantity,  pressure,  is  specified  while  the  others  are 

extrapolated  from  the  interior.    For  the  density  and  the  normal  velocity,  simple 

first-order  extrapolation  is  used,  and  the  axial  outflow  velocity  is  obtained  from 

the  zero  order  outgoing  Riemann  invariant.    The  outflow  boundary  conditions 

are  given  by 

ui  =  Ri"-2ai/(r-i).    ai=VWA 

Wj  =  w2 

Pl=p2 

(67) 
On  the  body  surface  the  nonslip  condition  is  applied  for  the  velocities.    The 
density  and  pressure  are  obtained  from  the  interior  by  extrapolation.   For  the  C- 
type  grids  used  in  this  study  averaging  of  the  flow  variable  at  the  wake  cut  is 
used. 

4.     Turbulence   Models 

a.     Baldwin-Lomax    Turbulence    Model 

The  Baldwin-Lomax  model  [Ref.  12]  is  a  two  -layer,  inner  and 
outer  eddy  viscosity  model  for  the  computation  of  two-  and  three-dimensional 
flows.  Patterned  after  the  Cebeci-Smith  model  [Ref.  11],  it  incorporates  a 
modification  that  bypasses  the  need  for  finding  the  edge  of  the  boundary  layer. 
This  is  achieved  by  introducing  the  vorticity  in  place  of  the  boundary  layer 
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thickness.    The  turbulent  effect  is  simulated  by  an  inner  eddy  viscosity,  given 
by: 


\r\  /inner       r^    r^l  J         ^crossover  //ro\ 

(68) 


where  ycrossover  is  the  smallest  value  of  y  at  which   the  values  of  the  inner  and 
outer  eddy  viscosities  are  the  same  and  the  Prandtl  mixing  length  1  is 

r  f-y^l 

l  =  kyil-exp  —   I 

L  VA    )\  (69) 

The  magnitude  of  the  vorticity  in  two-dimensions  Icol  and  y+  are  defined  as 
follows 


M= 


du    dv}  +_Pw"Ty  ^VPwuwY 


y  = 


dy     <?xj  and  M*  Mw  (jo) 

where  A+  is  a  constant.  The  subscript  w  denotes  values  at  the  wall  or  airfoil 
surface  in  this  case. 

The  outer  eddy  viscosity  is  given  by  the  following  expression 

vA)outer~^ClP^yJwAKE^y''raEB  y<:rossover  ~  Y  (1\\ 

where  k  and  Ccp  are  constants  and  F(y)wake=ymax  Fmax  for  boundary  layers  and 
F(y)wake=CWkymax  (UDIF2/Fmax  )  for  wakes  and  separated  boundary  layers. 
The  quantities  ymax  and  Fmax  are  determined  from  the  expression 


r  (-  +  YI 

F(y)  =  y  |o|    1-expl-^rj 


(72) 
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The  exponential  term  of  Equation  (72)  is  set  to  zero  for  wakes.  Fmax  is  the 
maximum  value  of  F(y)  that  occurs  in  a  profile  and  ymax  is  the  value  of  y  at 
which  it  occurs.  F(y)KLEB  is  the  Klebanoff  intermittency  factor  given  by 


F(y) 


r 


/, 


KLEB 


1  +  5.5 


^  y 


KLEBy 
man     / 


n 


The  quantity  Udif  is  the  difference  between  Umax=Uy=ymax 
total  velocity  in  the  profile 


(73) 
and  the  minimum 


UDIF=(vu^ 


+  v' 


y=y„: 


-W? 


+  v 


y=y, 


(74) 


The  second  term  in   Udif  is  taken  to  by  zero  except  in  the  case  of  wakes. 

The   constants   were   determined    to   achieve   agreement   with 
Reference    1 1    and  can   be  found   in   the  original  document.      This  model 
completely  models  turbulent  flow  over  an  airfoil,  it  is  relatively  easy  to  code,  it 
does  not  degrade  the  solution  convergence,  and  it  is  computationally  efficient. 
b.    Johnson-King   Turbulence   Model 

The  Johnson-King  model  [Ref.  13]  takes  into  account  the 
convective  and  diffusive  effects  on  the  Reynolds  shear  stress  -u'w'  in  the 
streamwise  direction.   The  eddy  viscosity  is  given  by 


r 


( ■ 


V.  =  V. 

1         lo 


1  -exp 


SI 


VVlo. 


(75) 


where  Vt ,  Vto  describe  the  eddy  viscosity  variation  in  the  inner  and  outer  part  of 
the    boundary    layer.        The    inner    eddy    viscosity    is    computed    as 
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vt  =  D2Ky^(-u'w') 
outer  viscosity  is  given  by 


2  ~(y/A+ ) 

max , where  D  =1  -e         ',  where  the  constant  A+=15.    The 


vto=a(x)[0. 168^5x7]  (76) 

where  7 is  the  Klebanoff  intermittency  function  7^[l+5.5(y/5)6]-l  and  g(x)  is  the 
solution  of  the  ordinary  differential  equation  which  describes  the  development  of 
-u'w'lmax  along  the  path  of  maximum  shear  stress.   This  model  accounts  for  the 
effects  of  convection  and  diffusion  on  the  Reynolds  stress  development  through 
the  solution  of  the  following  ODE 


ck      2lL,L, 


1- 


+ 


^dif^m 


'°q  J 


a4Q7-(f)J 


1- 


I  l      yon  j       (77) 

here  Cdif  and  ai  are  modeling  constants,  um  is  the  maximum  average  mean 
velocity  and  g=[-u,w'lmax]"1/2,  geq^t-u'w'lmax.eg]"172  where  Lm  is  the 
dissipation  length  evaluated  as 


Lm=Q40y  ym/5<Q225 

Lm=Q09y         ym/5>Q225 


(78) 
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The  equilibrium  shear  stress  geq  in  Equation  (77)  is  determined  from  the 
following  equilibrium  eddy  viscosity  distribution 


r 


V       =  V 

Ki,eq  %.eq 


1-exp 


ij.eq 
V 

V,  to.eq  J 


vtoeq  =<r(x)[Q168Uc5>]  (79) 

An  implicit  Euler  method  is  used  for  the  numerical  solution  of  Equation  (77), 
and  the  maximum  shear  stress  at  each  iteration  level  is  updated  as  follows 


n+1 
/     \n+l  /      \"    v> 

(7{X)  =<J{X) 


Vlo  (80) 


Solutions  with  the  Johnson-King  turbulence  model  were  obtained  as  follows. 
First  a  convergent  solution  using  the  Baldwin-Lomax  turbulence  model  for  the 
entire  flow  field  was  computed.  Then  the  Johnson-King  model  was  applied  only 
to  the  upper  part  of  the  airfoil.  To  initiate  the  solution  a(x)  in  Equation  (76)  is 
set  equal  to  unity  and  it  is  allowed  to  change  according  to  Equation  (80)  until 
the  final  solution  is  obtained. 

c.     Algebraic    RNG-based    Turbulence    Model 

Recently  an  algebraic  eddy  viscosity,  as  well  as  a  two-equation  k- 
£  model  based  on  Renormalization  Group  (RNG)  Theory  of  turbulence  [Ref.  14] 
were  proposed  for  the  closure  of  the  Reynolds-averaged  Navier-Stokes 
equations.  The  algebraic  model,  although  free  from  uncertainties  related  to  the 
experimental  determination  of  empirical  modeling  constants,  still  requires 
specifications  at  an  integral  length-scale  of  turbulence  which  reduces  the 
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generalities  of  the  model.  Here  the  integral  scale  is  assumed  proportional  to  the 
distance  from  the  wall  and  the  eddy  viscosity  for  the  RNG-based  algebraic 
turbulence  model  is  obtained  as  in  Martinelli  [Ref.  15]  from  the  following 
formula 


v  =  v, 


r     f»  (\     i  y4    n 


(81) 


where  v=Vt+vi,  H  is  the  Heavyside  step  function  and  <J)  is  the  dissipation 
function  (j)=Tij  (  du[  I  3xj  ).  The  RNG-based  turbulence  model  is  applied  only  for 
the  suction  surface  separated  flow  region  while  the  pressure  side  and  the  wake 
regions  are  computed  with  the  Baldwin-Lomax  model. 
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IV.     RESULTS    AND    DISCUSSION 

A.   METHOD   OF   INVESTIGATION 

As  stated  in  the  introduction  this  investigation  is  divided  into  three  parts; 
validation  of  the  Navier-Stokes  code  (Ns2.f)  through  comparison  with 
experimental  data  and  a  well  tested  inviscid  unsteady  Panel  method  code  as 
well  as  a  steady  viscous/inviscid  code.  Finally  evaluation  of  the  effects  of 
reduced  frequency,  mean  angle,  amplitude,  Mach  number,  and  Reynolds 
number  on  unsteady  flow  behavior  are  examined.  Also,  during  the  course  of 
this  study  it  became  necessary  to  investigate  the  effects  of  several  turbulence 
models  on  the  results  computed  for  the  harmonically  oscillating  cases. 
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The  following  table  lists  the  cases  studied  during  the  course  of  this 
investigation: 

TABLE   2.  TEST    CASES 


Case 

Type  of  Motion 

Parameters 

1 

Steady  State 
(Transonic) 

M=0.7,a=1.4°,Re=9xl06 

2 

Steady  State 
(Transonic) 

M=0.799,  a=2.26°,  Re=9xl()6 

3 

Steady  State 

M=0.3,Re=3xl()6 

4 

Steady  State 

M=0.3,Re=4xl06 

5 

Rapidly  Pitching  Airfoil 
(Ramp) 

M=0.3,  k=0.01272,  Re=2.7xl06 

6 

Harmonically  Oscillating 
(Multiple  Parameters) 

a(t)=Ao°±Ai°sin(cot) 

7 

Harmonically  Oscillating 
(Small  Amplitude) 

a(t)=13°±2.5°sin(G)t),  k=0.2,  Re=4xl()6 

8 

Harmonically  Oscillating 
(Small  Amplitude) 

a(t)=13°±2.5°sin(G)t),  k=0.1,  Re=2xl06 

9 

Harmonically  Oscillating 

a(t)=9°±5°sin(cot),  k=0.2,  Re=4xl06 

10 

Harmonically  Oscillating 

a(t)=10±5.5°sin(cot),  k=0.1,  Re=4xl06 

11 

Harmonically  Oscillating 

a(t)=ir±5°sin(a)t),  k=0.1,  Re=4xl06 
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1.     Steady-State 

Flow  solutions  can  be  obtained  either  by  rotating  the  grid  or  by  rotating 
the  flow.   Cases  1  and  2  solutions  were  obtained  by  rotating  the  grid. 

First  the  procedure  was  used  for  the  computation  of  steady  flow  over 
the  NACA  0012  airfoil.  Steady-state  solutions  were  obtained  for  subsonic  and 
transonic  flow  regime.  Flow  solutions  for  Case  1  and  2  (Table  2)  were 
computed  by  rotating  the  direction  of  the  flow  relative  to  the  stationary  grid. 
However,  in  Cases  3  and  4  the  grid  was  rotated  to  the  desired  angle  using  a 
simple  Fortran  code  called  rotgr.f  and  the  free  stream  flow  remained  parallel  to 
the  x-axis.  The  program  also  translates  the  grid  to  a  desired  pivot  point.  In  this 
thesis  the  pivot  point  was  always  chosen  to  be  the  quarter  chord  point.  Either 
method  produced  the  same  solution.  The  steady-state  solutions  tested  the 
accuracy  of  Ns2.f  and  provided  initial  solutions  for  the  unsteady  cases.  Ns2.f 
was  also  compared  to  an  existing  steady,  incompressible,  viscous/inviscid 
method  (Incompbl.f)  obtained  from  Dr.  T.  Cebeci,  California  State  University  at 
Long  Beach. 

a.     Case  1.    M=0.7,  a=1.4°,    Re=9xl06 

The  flow  conditions  for  the  first  test  case  are  Mx>=0.7,  at  an  angle 
of  attack  of  1.4  degrees,  and  a  Reynolds  number  of  nine  million.  Explicit 
boundary  conditions  were  used  and  the  Baldwin-Lomax  turbulence  model  was 
selected.  The  solution  was  obtained  on  a  161x64  point  grid.  This  case  was 
chosen  because  the  solution  converged  rapidly.  Figure  3  shows  good 
agreement  between  the  inviscid  (Euler)  pressure  distribution  and  the  measured 
values  of  Harris  [Ref.  16].  Next,  a  Navier-Stokes  solution  was  obtained  whose 
convergence  history  is  given  in  Figure  4.    The  solution  was  carried  out  to  2500 
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iteration  at  a  timestep  of  approximately  0.08.  The  pressure  distribution 
computed  by  Ns2.f  (Figure  5)  is  also  in  good  agreement  with  the  measured 
values.  However,  the  suction  peak  is  more  accurately  predicted  by  the  inviscid 
solution.  For  a  free  stream  Mach  number  (Moo)  of  0.7  a  maximum  Mach 
number  of  1.1  (Figure  6)  is  obtained  on  the  upper  surface  near  the  leading  edge. 
The  density  contour  plot,  Figure  7,  shows  the  smooth  contours  indicative  of  a 
fully  converged  solution.  The  velocity  field,  Figure  8,  shows  the  boundary  layer 
thickness,  represented  by  the  thick  dark  line.  The  sonic  line  is  represented  by 
the  line  extending  from  7%  chord  to  20%  chord. 
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Figure  3.     Pressure  Distribution  (Inviscid),  M=0.7,  a=1.4°,   Re=9xl06 


39 


* 

O 

X 

On 

II 

<D 

#^ 

d 


>* 


^ 


o 

c 

b/j 
> 

a 

o 

U 


o 
o 
o 
m 


o 
o 
o 

(N 


o 
o 
o 


r~~ 


oo 


On 


O 

o 


c 

o 

a 

5- 


jnnpis^y  xuj/\[ 


Figure  4.     Convergence  History,  M=0.7,  <x=1.4\   Re=9xl06 
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Figure  5.     Pressure  Distribution  ,  M=0.7,  a  =  1.4°,   Re=9xl06 
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Figure  6.    Mach  Contour,  M=0.7,  a=1.4°,   Re=3xl0<> 
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Figure  7.     Density  Contour,  M=0.7,  a=1.4°,   Re=9xl06 
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Figure  8.    Velocity  Field,  M=0.7,  a=1.4°,   Re=9xl06 
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b.     Case  2.     M=0.7  99,  a=2.26°,    Re=9xl0^ 

The  next  steady-state  case  investigated  was  at  a  Mach  number 
equal  to  0.799,  at  an  angle  of  attack  of  2.26  degrees  and  a  Reynolds  number  of 
nine  million.  As  in  Case  1,  an  Euler  solution  was  first  obtained  and  compared  to 
measured  data  as  well  as  to  a  viscous  solution  obtained  with  an  upwind  scheme 
using  the  Johnson-King  turbulence  model.  As  seen  in  Figure  9,  the  pressure 
distribution  for  the  inviscid  solution  is  in  fair  agreement  with  the  measured  data 
until  the  shock  is  reached  at  the  50%  chord  point.  After  this  point  the  computed 
solution  fails  to  predict  the  pressure  jump  due  to  the  presence  of  a  shock  at 
mid-chord.  Note  that  the  upwind  scheme  in  combination  with  the  Johnson-King 
turbulence  model  accurately  predicts  the  location  of  the  shock,  even  though  the 
magnitude  of  the  pressure  gradient  is  off  slightly.  Next,  a  viscous  solution  was 
computed  by  Ns2.f  with  the  Baldwin-Lomax  turbulence  model.  Fourth  order 
dissipation  was  added  to  produce  a  smooth  flow  solution  across  the  weak  shock 
located  at  mid  chord.  The  solution  was  continued  for  6000  iterations.  The 
spike  in  the  convergence  history  (Figure  10)  was  due  to  the  addition  of  more 
fourth  order  dissipation.  It  is  seen  that  the  shock  location  was  not  predicted. 
At  this  angle  of  attack  and  Mach  number  the  flow  starts  to  separate  and  the 
Baldwin-Lomax  turbulence  model  fails  to  model  this  separated  flow  region. 
Figure  1 1  shows  the  viscous  pressure  distribution  comparison.  The  sonic  line  is 
shown  on  the  Mach  contour  plot  of  Figure  12.  A  maximum  Mach  Number  of 
1.45  is  reached  in  this  supersonic  flow  region.  Note  the  steep  Mach  gradient 
(Figure  12)  at  approximately  60%  chord.  Figures  13  and  14  show  the  pressure 
contour  and  density  contour,  respectively.  Again  the  steep  density  gradient 
indicates  the  shock  location  at  approximately  60%  chord  vice  the  measured 
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location  at  50%  chord.    Finally,  the  velocity  field  and  sonic  line  are  shown  in 
Figure  15. 
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Figure  9. 


Pressure  Distribution  (Inviscid)  ,  M=0.799,  a  =  2.26°, 
Re=9xl06 
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Figure  10.     Convergence  History  (Ns2.f),  M=0.799,  a=2.26°, 

Re=9xl06 
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Figure  11.     Pressure  Distribution,  M=0.799,  a=2.26°,   Re=9xl06 


49 


-   -   =  -    - 


s  3  s  =  =  =   .=  3333333333  = 

iiiiiZ  -:5§§§8888§§  = 


-33 N 

Z 

3  •  3  3  3  3 


Figure  12.     Mach  Contour,  M=0.799,  a=2.26°,    Re=9xl()6 
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Figure  13.     Pressure  Contour,  M=0.799,  a=2.26°,    Re=9xlO<5 
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Figure  14.     Density  Contour,  M=0.799,  a=2.26°,    Re=9xl()6 
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c.     Case  3.     M=0.3,  Re=3xl06 

In  order  to  produce  C\  vs  a  and  Cm  vs  a  curves,  solutions  were 
computed  for  a=2°,  a=4°,  a=6°,  a=8°,  a=9°,  a=10°,  a=12°,  a=13°,  a=13.5°, 
a=14°,  and  a=15°.    However,  only  results  for  a=4°,  a=ll°,  and  ct=14°,  are 
presented. 

Figure  16  shows  the  convergence  history  of  Ns2.f  for  a  161x64 
algebraic  C-type  grid.  The  Courant-Friedrichs-Levy  (CFL)  condition  was 
approximately  2000-2500  which  corresponded  to  a  timestep  between  0.005  and 
0.01.  A  comparison  of  the  computed  pressure  distributions  at  a=4°  (Figure  17) 
shows  good  agreement  between  the  Navier-Stokes  calculations  using  the 
Baldwin-Lomax  turbulence  model  and  the  viscous/inviscid  code  predictions. 
The  measured  suction  peak  is  somewhat  higher  than  predicted  by  the  codes. 
Figures  18  and  19  show  fair  skin  friction  agreement.  Note  that  the  Navier- 
Stokes  code  assumes  turbulent  flow  over  the  entire  airfoil,  while  the 
viscous/inviscid  interaction  code  uses  Michel's  method  to  predict  the 
laminar/turbulent  transition.  As  seen  in  Figure  18,  the  transition  point  for  the 
upper  surface  is  at  10%  chord.  A  Mach  contour  plot  is  shown  in  Figure  20  and 
a  density  contour  plot  with  the  normalized  stagnation  pressure  (0.98)  contour 
overlayed  in  Figure  21.  A  maximum  Mach  number  of  0.48  was  attained.  The 
density  contours  are  smooth  with  steep  gradients  at  the  leading  edge. 

In  Figures  22  through  29  similar  results  are  presented  for  angles  of 
attack  of  11  and  14  degrees.  It  is  seen  that  the  agreement  between  the  two 
codes  and  the  measured  pressure  distributions  is  quite  satisfactory,  although  the 
suction  peak  tends  to  be  underpredicted  by  both  codes. 


54 


Figures  30  and  31    show  the  Civs  a  and  Cm  vs  a    curves, 
respectively.    Stall  is  seen  to  occur  at  approximately  a=13.5°.    The  Navier- 
Stokes  code  reproduces  the  experimental  lift  and  moment  values  up  to  the 
maximum  lift  value  quite  well,  whereas  the  viscous/inviscid  code  shows  greater 
deviation  for  angles  of  attack  exceeding  about  eight  degrees. 
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Figure  16.     Convergence  History,  M=0.3,  a=4°,   Re=3xl06 
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Figure  17.     Pressure  Distribution,  M=0.3,  a=4°,   Re=3xl06 
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Figure  18.     Skin  Friction  Distribution  (Upper  Surface),  M=0.3,  a=4°, 

Re=3xl06 
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Figure  19.     Skin  Friction  Distribution  (Lower  Surface), 

Re=3xl06 


M=0.3,  a  =  4°, 


59 


c 

X  _  = 

a.  r  - 

_  t    -  a 

x  at    *-  a 


S     IB 


atc3a  =  =  r=  =  =  =  =  =  =r       rxxcoexcxx 

Z     .    .  "■  ~  ~  ...  ~  '.".'".     '.  I'.'..'.'.'.'.. 

C^—    ~    —    ------    ~    —    ~-  r    —    =S3C3XX    = 


Figure  20.     Mach  Contour,  M=0.3,  a=4°,   Re=3xl06 
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Figure  21.     Density  Contour,  M=0.3,  a=4°,   Re=3xl06 
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Figure  22.     Convergence  History,  M=0.3,  a=ll°,    Re=3xl06 
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Figure  23.     Pressure  Distribution,  M=0.3,  a=ll°,   Re=3xl06 
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Figure  24.       Skin  Friction  Distribution  (Upper  Surface),  M=0.3, 

a=ll°,   Re=3xl06 
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Figure  25.       Skin  Friction  Distribution  (Lower  Surface),  M=0.3, 

a=ll°,   Re=3xl06 
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Figure  26.     Mach  Contour,  M=0.3,  a=ll°,   Re=3xlQ6 
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Figure  27.     Density  Contour,  M=0.3,  a=ll°,   Re=3xl06 
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Figure  28.     Convergence  History,  M=0.3,  a=14°,    Re=3xl06 
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Figure  29.     Pressure  Distribution,  M=0.3,  a=14°,   Re=3xl0<> 
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Figure  30.    C|  vs  a,  M=0.3,  Re=3xl0<> 
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Figure  31.    Cm  vs  a  ,  M=0.3,  Re=3xl()6 
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d.     Case  4.     M=0.3,  Re=4xl06 

Figure  32  shows  a  comparison  of  the  computed  lift  curve  with  the 
experimental  data  at  a  Reynolds  number  of  four  million.  It  can  be  seen  that  the 
computed  maximum  lift  angle  and  hence  the  static  stall  angle  is  13.5°  degrees. 
This  agrees  quite  well  with  the  experimentally  measured  stall  angle  of 
McCroskey  et  al  [Ref.  17].  The  computed  lift  and  pitching  moment  are  in  fair 
agreement  for  small  angles  of  attack,  but  start  to  deviate  at  larger  angles. 
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M  =  0.3,  Re  =  4  x  106  Baldwin-Lomax  Turbulence  Model 

Experiment:  McCroskey 
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Figure  32.    Ci  vs.  a  and  Cm  vs.  a,  M=0.3,  Re=4xl06 
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2.      Rapidly  Pitching  Airfoil 

Next,  AGARD  Case  6  from  Landon  [Ref.  18]  was  investigated  using 
Ns2.f  and  U2diif.f.  The  computed  results  are  compared  to  the  measured  data. 
The  flow  condition  for  this  case  corresponds  to  a  freestream  Mach  number 
equal  to  0.3  at  a  non-dimensional  pitch  rate  equal  to  0.01272  and  a  Reynolds 
number  of  2.7  million.  Flow  solutions  for  a  rapidly  pitching  airfoil  were 
computed  by  pitching  the  airfoil  about  the  quarter  chord  at  a  constant  rate  from 
zero  degrees  angle  of  attack  to  a  final  angle  of  attack  of  15.54°.  Figure  33 
shows  a  sketch  of  the  motion  produced  by  a  rapidly  pitching  airfoil. 
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Figure  33.     Rapidly  Pitching  Airfoil  Motion 

First  a  steady-state  solution  at  zero  degrees  angle  of  attack  was  well 
converged  to  3800  iterations.    Then  the  iteration  counter  was  reset  to  zero  and 
the  airfoil  was  rapidly  pitched  up  to  the  final  angle.   The  reduced  frequency  k  is 
given  by  k=oc  c/  2  Uoo,  where  a  is  the  pitch  rate.    In  terms  of  non-dimensional 
quantities  a(t)=2Mook  and  cx(t)=  cto+ai(t/  T),  where  ao  (Aq),  is  the  initial  angle 
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of  attack  a=0°,  a\  (Aj),  is  the  final  angle  of  attack  a=15.54°,  and  T  is  the  time 
required  for  the  airfoil  to  complete  the  pitching  motion. 

a.     Case  5.     M=0J,     k=0.01272,  Re=2Jxl()6 

Figure  34  shows  the  comparison  of  the  computed  pressure 
distribution  for  the  ramp-type  unsteady  motion  at  a=4.84°,  for  an  inviscid  Euler 
calculation.  This  Euler  solution  uses  a  121x35  C-grid  for  the  computation. 
Similar  results  were  produced  for  ensuing  angles  of  attack.  Figures  35  through 
39  show  pressure  distributions  at  a=4.84°,  a=6.72°,  a=10.80°,  a=12.83°,  and 
oc=15.54°.  Surprisingly,  the  inviscid  code  predicted  the  suction  peak  more 
accurately  than  the  N-S  code  for  all  angles  of  attack.  Several  explanations  are 
offered  to  account  for  this  unexpected  result.  First  the  161x65  C-grid  may  not 
be  fine  enough  to  predict  the  magnitude  of  the  suction  peak  and  the  N-S  code 
may  have  too  much  numerical  viscosity,  thus  in  effect  performing  calculation  at 
a  lower  Reynolds  number.  The  flow  at  the  leading  edge  may  in  fact  be  laminar, 
while  the  computed  solution  assumed  turbulent  flow  everywhere.  Figure  40 
shows  the  Mach  contour  and  Figure  41  shows  the  density  contour  at  oc=10.8°. 
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Figure  34.     Pressure  Distribution  (Inviscid),  Ramp  Motion,  M=0.3, 

<x=4.87°,   k=0.01272,   Re=2.7xl()6 
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Figure  35.     Pressure  Distribution,  Ramp  Motion,  M=0.3,  oc=4.87 

k=0.01272,    Re=2.7xl06 
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Figure  36.     Pressure  Distribution,  Ramp  Motion,  M=0.3,  ct=6.72' 

k=0. 01272,    Re=2.7xl0<> 
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Figure  37.     Pressure  Distribution,  Ramp  Motion,  M=0.3,  a  =  10.80°, 

k=0.01272,    Re=2.7xl06 
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Figure  38.     Pressure  Distribution,  Ramp  Motion,  M=0.3,  a=  12.83°, 

k=0.01272,    Re=2.7xl06 
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Figure  39.     Pressure  Distribution,  Ramp  Motion,  M=0.3,  a=15.54 

k=0.01272,    Re=2.7xl06 
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Figure  40.    Mach  Contour,  Ramp  Motion,  M=0.3,  a=10.8°, 

k=0.01272,    Re=2.7xl0* 
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Figure  41.     Density  Contour,  Ramp  Motion,  M=0.3,  cc=10.8°, 

k=0.01272,    Re=2.7xl06 
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A  detail  of  the  flow  field  generated  by  Ns2.f  at  a=  12.83°  is  shown 
in  Figure  42.  The  flow  appears  to  be  smooth  and  attached  (the  dark  line  shows 
the  boundary  layer  edge).  However,  magnification,  as  shown  in  Figure  43, 
reveals  the  start  of  a  slight  reverse  flow  region  at  approximately  10%  chord.  At 
a=15.54°,  Figure  44,  the  boundary  layer  is  much  thicker  with  reverse  flow 
profiles  from  approximately  10%  chord  to  the  trailing  edge. 

Figures  45  and  46  present  the  computed  C\  vs  a  and  Cm  vs  a 
curves,  respectively.  The  Navier-Stokes  code  underpredicts  lift  and  pitching 
moment  coefficients  whereas  the  inviscid  panel  method  gives  significantly 
better  lift  predictions.  However,  it  fails  to  predict  the  pitching  moment 
coefficient  for  angles  greater  than  10°.  A  more  comprehensive  investigation  of 
the  dynamics  of  a  rapidly  pitching  airfoil  can  be  found  in  Grohsmeyer  [Ref.  19]. 
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Figure  42.     Velocity  Field,  Ramp  Motion,  M=0.3,  a=12.83°, 

k=0.01272,    Re=2.7xl06 
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Figure  43.     Velocity  Field  (Magnified),  Ramp  Motion,  M=0.3, 
a=12.83°,   k=0.01272,   Re=2.7xl0<> 
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Figure  44.    Velocity  Field,  Ramp  Motion,  M=0.3,  a=15.54 

k=0. 01272,    Re=2. 7x106 
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Figure  45.    Ci  vs  a,  Ramp  Motion,  M=0.3,  k=0.01272,  Re=2.7xl06 
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Figure  46.    Cm  vs  a  ,  Ramp  Motion,  M=0.3,  k=0.01272,  Re=2.7xl06 
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3.      Harmonic  Airfoil  Oscillation  Using  the  Baldwin-Lomax 
Turbulence   Model 

Next,  the  code  was  applied  to  sinusoidal  airfoil  oscillations  about  the 
quarter  chord  point.  These  time  periodic  solutions  were  obtained  for  the  second 
cycle  because  the  results  obtained  for  the  third  cycle  were  indistinguishable 
from  ones  of  the  second  cycle.  Figure  47  shows  the  computed  angle  of  attack 
history. 
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Figure  47.     Harmonically  Oscillating  Airfoil 


Here,  Ao  is  the  mean  angle  of  attack,  and  A]  is  the  amplitude.  For 
each  oscillatory  case  a  steady-state  solution  was  obtained  for  the  minimum 
angle  of  attack  during  the  cycle.  The  iteration  counter  was  set  to  zero  as  the 
sinusoidal  motion  was  time  shifted  half  a  cycle  so  that  the  start  of  the  motion 
begins  at  the  minimum  angle  of  attack.  This  ensures  a  smooth  transition  for  the 
converged  steady  state  solution  to  the  start  of  the  oscillatory  motion.     Ten 
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thousand  (10,000)  iterations  were  computed  per  cycle,  with  a  timestep  equal  to 
approximately  0.05  to  0.01. 

a.    Case  6.    a(t)=Ao°±Aj°sin((Ot) 

All  the  computations  presented  in  this  section  are  computed  with 
Ns2.f  using  the  Baldwin-Lomax  turbulence  model.  Figures  48  and  49  show  the 
effect  of  Mach  number  on  the  lift  and  moment  hysteresis  loops  for  small 
amplitude  oscillations  about  a  mean  angle  of  12  degrees  with  an  amplitude  of 
3.0  degrees,  at  a  Reynolds  number  of  four  million  and  at  a  reduced  frequency 
of  0.1.  Note  that  the  hysteresis  loops  for  Moo=0.4  oscillate  wildly  and  third 
cycle  computations  do  not  coincide  with  results  from  previous  cycles.  For  the 
Moo=0.3  computation  the  moment  hysteresis  loop  indicates  a  stable  condition. 
The  effect  of  amplitude  is  shown  in  Figures  50  and  51.  Figures  52  and  53  also 
show  the  effect  of  amplitude  for  a  larger  mean  angle  of  attack,  Ao=14.0\  at  the 
same  Mach  number,  reduced  frequency,  and  Reynolds  number,  as  before.  Both 
solutions  oscillate  in  a  manner  that  is  not  repeatable  if  the  motion  is  allowed  to 
continue  for  additional  cycles.  For  this  comparison  the  maximum  angle  of 
attack  attained  during  a  cycle  exceeded  16.5°.  A  comparison  of  the  first  and 
second  cycle  lift  and  moment  coefficient  loops  is  shown  in  Figures  54  and  55. 
Note  that  the  second  cycle  lift  curve  starts  at  a  slightly  greater  lift  coefficient 
than  in  the  first  cycle.  Aside  from  this  anomaly,  the  loops  are  in  good 
agreement  but  exhibit  no  instability. 

Figure  56  shows  the  computed  lift  and  moment  loops  for  small 
amplitude  oscillations  about  a  mean  angle  of  13  degrees  with  an  amplitude  of 
2.5  degrees,  at  a  Reynolds  number  of  two  million,  for  three  values  of  reduced 
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frequency.    It  can  be  seen  that  all  the  moment  loops  computed  in  this  section 
using  the  Baldwin-Lomax  turbulence  model  indicate  torsional  stability. 
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Figure  48.     Effect  of  Mach  Number  on  Unsteady  Motion, 
a(t)=12°±3°sin(cot),  k=0.10,  Re=4xl0<> 
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Figure  49.     Effect  of  Mach  Number  on  Unsteady  Motion, 
a(t)=12°±3°sin(cot),  k=0.10,  Re=4xl()6 
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Figure  50.     Effect  of  Amplitude  on  Unsteady  Motion,  M=0.3, 
a(t)=12o±Ai°sin((0t),  k=0.10,  Re=4xlQ6 
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Figure  51.     Effect  of  Amplitude  on  Unsteady  Motion,  M=0.3, 
a(t)=12°±Ai°sin(G>t),  k=0.10,  Re=4xl06 
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Figure  52.     Effect  of  Amplitude  on  Unsteady  Motion,  M=0.3, 
a(t)  =  14o±Aiosin(0)t),  k=0.10,  Re=4xl0<> 
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Figure  53.     Effect  of  Amplitude  on  Unsteady  Motion,  M=0.3, 
a(t)  =  14°±Ai°sin(cot),  k=0.10,  Re=4xl0<> 
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Figure  54.    Effect  of  1st  and  2nd   Cycle  on  Unsteady  Motion,  M=0.3, 
a(t)=13°±  2°sin(cot)/k=0.10,   Re=4xl()6 
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Figure  55.     Effect  of  1st  and  2nd   Cycle  on  Unsteady   Motion,  M=0.3, 
a(t)  =  13±2°sin(cot),  k=0.10,  Re=4xl06 
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Figure  56.         Effect  of  Reduced  Frequency  on  Unsteady  Motion 
M=0.3,  a(t)  =  13±2.5°sin(cot),   Re=2xl<)6 
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b.     Case  7.    a(t)=130±2.5°sin((Dt),   k=0.2,  Re=4xlO^ 

In  an  effort  to  determine  if  the  flow  field  was  in  fact  predicting  a 
trailing  edge  vortex  that  would  shed  during  the  oscillation  cycle  as  seen  in  the 
experimental  results,  instantaneous  particle  traces  and  the  corresponding 
velocity  fields  were  plotted  for  the  downstroke  at  a  Mach  number  of  0.3,  a 
mean  angle  of  attack  of  13  degrees,  an  amplitude  of  2.5  degrees,  a  reduced 
frequency  of  0.1  and  a  Reynolds  number  of  four  million.  Figures  57  through  61 
show  the  downstroke  portion  of  the  cycle  beginning  at  a=15.3°  and  ending  at 
a=14.5°.  A  recirculatory  region  is  seen  in  Figure  57  which  grows  slightly,  as 
seen  in  Figure  58.  This  region  then  appears  to  diminish  in  size  and  intensity. 
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Figure  57.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  a=15.3°   (Downstroke), 

a(t)=13°±2.5°sin(cot),  k=0.20,   Re=4xl06 
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Figure  58.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  <x=15.1°   (Downstroke), 

a(t)=13°±2.5°sin(cot),  k=0.20,  Re=4xl06 
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Figure  59.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  a=14.9°    (Downstroke), 

a(t)=13°±2.5°sin((0t),  k=0.20,  Re=4xl06 
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Figure  60.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  a=14.7°   (Downstroke), 

a(t)  =  13°±2.5°sin(cot),  k=0.20,  Re=4xl06 
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Figure  61.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  a=14.5°   (Downstroke), 

a(t)  =  13°±2.5°sin(G)t),  k=0.20,  Re=4xl06 
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c.    Case  8.    a(t)=13°±2.5°sin(cot),   k=0.1,  Re=2xl0^ 

The  plots  generated  for  Case  7  did  show  a  recirculatory  region, 
but  a  full  understanding  of  its  development  and  structure  could  not  be  gained 
from  only  a  portion  of  the  oscillatory  cycle.  Therefore  in  Figures  62  through  80 
a  complete  cycle  is  shown  for  the  same  flow  condition  as  in  Case  7  with  only 
the  reduced  frequency  being  increased  from  0.1  to  0.2  and  Reynolds  Number 
decreased  from  four  million  to  two  million.  A  recirculatory  region  is  seen  to 
form  at  a=14.9°  on  the  upstroke,  which  continues  to  grow  and  intensify.  This 
recirculatory  region  reaches  its  largest  value  at  cc=15.3°  on  the  downstroke  and 
then  diminishes  and  finally  the  flow  becomes  completely  reattached.  It  is 
important  to  note  that  even  for  this  relatively  fast  oscillation  the  recirculatory 
region  did  not  shed  into  the  wake. 
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Figure  62.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  a=12°  (Upstroke),  a(t)  =  13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  63.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  a=13°  (Upstroke),  a(t)=13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  64.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  ct=13.5°  (Upstroke),  a(t)  =  13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  65.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  a=14°  (Upstroke),  a(t)=13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  66.     Velocity  Field,  Oscillatory  Motion,  M=0.3,  a=14.5 
(Upstroke),   a(t)  =  13°±2.5°sin(a)t),  k=0.10,   Re=2xl0<> 
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Figure  67.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  a=14.7°  (Upstroke),  a(t)  =  13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  68.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  a=14.9°  (Upstroke),  a(t)  =  13°±2.5°sin((Ot)! 

k=0.10,    Re=2xl06 
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Figure  69.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  <x=15.1°  (Upstroke),  a(t)  =  13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  70.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  a=15.3°  (Upstroke),  a(t)  =  13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  71.     Instantaneous  Particle  Trace  and  Velocity  Field, 
Oscillatory  Motion,  M=0.3,  <x=15.5°  (Upstroke),  a(t)=13°±2.5°sin(cot), 

k=0.10,    Re=2xl06 
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Figure  72.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  <x=15.3°   (Downstroke), 

a(t)  =  13°±2.5°sin(cot),  k=0.10,  Re=2xlQ6 
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Figure  73.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  a=15.1°   (Downstroke), 

a(t)=13°±2.5°sin(cot),  k=0.10,  Re=2xl06 
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Figure  74.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  ot=14.9°   (Downstroke), 

a(t)=13°±2.5°sin(cot),  k=0.10,  Re=2xl0<> 
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Figure  75.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  <x=14.7°   (Downstroke), 

a(t)=13°±2.5°sin(Q)t),  k=0.10,  Re=2xl06 
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Figure  76.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  a=14.5°   (Downstroke), 

a(t)=13°±2.5°sin(cot),  k=0.10,  Re=2xl06 
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Figure  77.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  a=14°    (Downstroke), 

a(t)  =  13°±2.5°sin(G)t),  k=0.10,  Re=2xl06 
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Figure  78.     Velocity  Field,  Oscillatory  Motion,  M=0.3,  a=13.5 
(Downstroke),   a(t)  =  130±2.5°sin(cot),  k=0.10,  Re=2xl0<> 
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Figure  79.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  <x=13°   (Downstroke), 

a(t)=13°±2.5°sin(cot),  k=0.10,  Re=2xl0<> 
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Figure  80.     Instantaneous  Particle  Trace  and  Velocity  Field, 

Oscillatory  Motion,  M=0.3,  <x=12°   (Downstroke), 

a(t)=13°±2.50sin(cot),  k=0.10,  Re=2xl06 
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4.      Effect  of  Turbulence  Modeling 

a.    Case  9.    a(t)=9°±5°sin((Ot),   k=0.2,  Re=4xl0^ 

Figures  81  and  82  show  a  comparison  of  the  lift  and  moment 
hysteresis  loops  using  the  Baldwin-Lomax  turbulence  model  for  an  oscillation  of 
5  degrees  amplitude  about  a  mean  angle  of  attack  of  9.0  degrees  with 
McCroskey's  experimental  data  for  a  Mach  number  of  0.3,  a  Reynolds  number 
of  four  million,  and  a  reduced  frequency  k=0.2.  It  is  seen  that  there  is  a 
substantial  difference  between  the  computed  and  the  experimental  lift 
hysteresis  loops.  This  difference  becomes  even  more  pronounced  when  the 
computed  and  experimental  moment  hysteresis  loops  are  compared  with  each 
other.  In  an  effort  to  understand  the  failure  of  Ns2.f  to  predict  the  experimental 
pitching  moment  hysteresis  loop,  the  pressure  distributions  for  several  upstroke 
and  downstroke  angles  of  attack  were  compared  to  measured  data  from 
Reference  1.  Figures  83  through  87  show  pressure  distributions  for  oc=11.9°  on 
the  upstroke  portion  of  the  cycle;  and  a=13.7°,  a=13.0°,  a=10.5°,  and  oc=8.9°  on 
the  down  stroke  cycle.  U2diif.f  overpredicts  the  suction  peak  at  oc=11.9°  on  the 
upstroke  and  Ns2.f  underpredicts  the  peak.  Even  though  the  values  are  off,  the 
general  shape  of  the  distribution  is  well  predicted.  However,  on  the  downstroke 
at  a=13.7°,  a  significant  difference  between  the  computed  and  the  measured 
distributions  is  observed.  The  measured  pressure  distribution  shows  a  lower 
suction  near  the  leading  edge  but  higher  suction  near  the  trailing  edge.  Neither 
of  these  changes  are  predicted  by  Ns2.f.  As  the  oscillation  continues  on  the 
downstroke  the  measured  distribution  tends  to  flatten  out  creating  a  "plateau". 
Just  after  a=8.9°  on  the  downstroke,  the  flow  begins  to  reattach  and  the  positive 
pressure  difference  at  the  trailing  edge  becomes  smaller.   The  failure  of  Ns2.f  to 
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predict  the  experimental  loop  results  from  the  inability  of  Ns2.f  to  accurately 
compute  the  pressure  distributions.  Since  the  pitching  moment  is  an  integrated 
quantity,  the  size  of  the  area  under  the  pressure  distribution  aft  of  the  quarter 
chord  must  be  greater  than  the  area  forward  of  the  quarter  chord  to  compute  a 
negative  or  downward  pitching  moment.  Ns2.f  clearly  does  not  predict  the 
experimental  pressure  distributions  during  the  downstroke. 
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Figure  81.    Ci  vs  a  ,  Oscillatory  Motion,  M=0.3,  a(t)=9°±5°sin(cot), 

k=0.2,    Re=4xl0<> 
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Figure  82.    Cm  vs  a  ,  Oscillatory  Motion,  M=0.3,  a(t)=9°±5°sin(cot), 

k=0.2,   Re=4xl0<> 
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Figure  83.     Pressure  Distribution,  Oscillatory  Motion,  M=0.3, 
a=11.9°  (Upstroke),  a(t)=9°±5°sin(o)t),  k=0.2  Re=4xl()6 


132 


o 

X 

II 

• 

II 

J* 

in 
II 

< 

#^ 

II 

O 

< 

c- 

m 

• 

o 
II 


a 
o 


o 

o 


fc£    C  .       ^ 
o      >-       .  ■  — 

—  -p.  -o  -a 

I    =:  3  3 

«  £-2- 

-  n  E  5 

c  o 


00 


^O 


^T 


r) 


00 


O 


O 


o 


o 
d 


c 

_o 

o 

o 


T3 
O 

U 


]U9I01JJ9CQ  3J11SS3JJ 


Figure  84.     Pressure  Distribution,  Oscillatory  Motion,  M=0.3, 
a=13.7°   (Downstroke),  a(t)=9°±5°sin(a)t),  k=0.2,  Re=4xl()6 
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Figure  85.     Pressure  Distribution,  Oscillatory  Motion,  M=0.3, 
a=13.0°   (Downstroke),  a(t)=9°±5°sin(cot),  k=0.2,  Re=4xl06 
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Figure  86.     Pressure  Distribution,  Oscillatory  Motion,  M=0.3, 
a=10.5°   (Downstroke),  a(t)=9°±5°sin(cot),  k=0.2,  Re=4xl0* 
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Figure  87.     Pressure  Distribution,  Oscillatory  Motion,  M=0.3,  a=8.9 
(Downstroke),   a(t)=9°±5°sin(cot),  k=0.2,  Re=4xl()6 
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b.     Case  10.    a(t)=10°±5.5°sin((Ot),   k=0.1,  Re=4xl()6 

The  measured  loop  consists  of  two  subloops,  a  clockwise 
"unstable"  subloop  and  a  counterclockwise  "stable"  subloop.  It  is  seen  that  the 
Navier-Stokes  calculations  in  combination  with  the  Baldwin-Lomax  turbulence 
model  fails  to  predict  the  destabilizing  clockwise  pitching  moment  loop.  For  this 
reason  the  sensitivity  of  the  computed  loops  to  different  turbulence  models  was 
studied  next.  Figures  88  and  89  show  the  computed  lift  and  pitching  moment 
loops  for  oscillation  about  a  slightly  higher  mean  angle  of  10  degrees  at  an 
amplitude  of  5.5  degrees.  The  Mach  and  Reynolds  numbers  are  again  0.3  and 
four  million,  and  k  was  decreased  to  0.1  to  minimize  the  effects  of  reduced 
frequency.  It  can  be  seen  that  the  Baldwin-Lomax,  Johnson-King,  and  the 
RNG  turbulence  models  produce  significantly  different  hysteresis  loops.  The 
RNG  model  produces  relatively  good  agreement  with  the  measured  lift 
hysteresis  loop  but  fails  again  to  predict  the  destabilizing  moment  subloop. 
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Figure  88.    C|  vs  a  ,  Oscillatory  Motion,  M=0.3, 
a(t)=10°±5.5°sin((Ot),  k=0.10,  Re=4xl()6 
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Figure  89.    Cm  vs  a  ,  Oscillatory  Motion,  M=0.3, 
a(t)=10°±5.5°sin(cot),  k=0.10,  Re=4xl()6 
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c.     Case  11.    a(t)=llo±5°sin((0t),   k=0.1,  Re=4xl0^ 

A  solution  was  also  computed  for  oscillations  about  a  mean  angle 
of  1 1  degrees  at  an  amplitude  of  5  degrees,  the  other  parameters  left 
unchanged.  Figures  90  and  91  show  the  computed  lift  and  moment  hysteresis 
loops.  As  before,  the  calculations  fail  to  predict  the  expected  destabilizing 
moment  subloop.  The  reason  for  this  failure  can  be  better  appreciated  by 
visualizing  the  separated  flow  structures  computed  with  these  turbulence 
models.  Figures  92  and  93  show  the  computed  flowfield  with  particle  traces 
and  the  velocity  field.  Figure  92  corresponds  to  the  instant  when  the  airfoil 
oscillates  through  an  incidence  of  15  degrees  during  the  downstroke.  Different 
turbulence  models  produce  substantially  different  recirculatory  flow  patterns  on 
the  upper  surface  near  the  trailing  edge.  The  Baldwin-Lomax  model  produced 
the  smallest  recirculatory  region  while  the  RNG  model  produced  the  largest.  In 
addition,  Figure  93  (a=15°)  shows  the  relative  magnitudes  of  the  velocity 
vectors.  Note,  near  the  surface  at  the  trailing  edge,  the  magnitude  of  the 
velocity  vectors  in  the  reverse  flow  region  calculated  by  the  Baldwin-Lomax 
model  are  much  smaller  than  the  region  calculated  by  the  Johnson-King  model 
or  RNG  model.  As  the  cycle  continues,  Figures  94  and  95  (ot=146)  show  the 
recirculation  region  to  shrink  slightly  and  loose  intensity  for  the  Baldwin-Lomax 
model,  while  the  regions  as  calculated  by  the  Johnson-King  and  RNG  continue 
to  expand  and  intensify.  As  noted  in  Case  8  the  recirculatory  flow  region 
grows  and  then  vanishes  again  during  the  course  of  the  oscillation.  Even 
though  the  structure  at  the  trailing  edge  resembles  a  vortex  it  does  not  shed 
from  the  trailing  edge.  The  occurrence  of  the  destabilizing  moment  loop 
appears  to  be  closely  associated  with  the  vortex  shedding  from  the  trailing 
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edge,  an  event  which  occurs  soon  after  the  static  angle  is  substantially 
exceeded  during  part  of  the  cycle.  Although  the  static  stall  angle  is  significantly 
exceeded  in  the  present  calculations  the  numerical  procedure  along  with  the 
turbulence  modeling  used  in  the  calculations  appears  to  be  unable  to  produce 
the  anticipated  shedding  process. 
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Figure  90.    Cj  vs  a 


Oscillatory  Motion,  M=0.3,  a(t)  =  ll°±5°sin(Q)t), 
k=0.10,    Re=4xl06 
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Figure  91.    Cm  vs  a  ,  Oscillatory  Motion,  M=0.3,  a(t)  =  ll°±5°sin(cot), 

k=0.10,    Re=4xl06 
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EFFECT  OF  TURBULENCE  MODELING  ON  THE  COMPUTED 
TRAILING  EDGE  FLOWFIELD 

M^  =  0.3,  a  (t)  =  11°  +  5~  sin  (cot),  k  =  0.1,  Re  =  4  x  106 

a  =  15    DOWNSTROKE 


(a)  Baldwin-Lomax 


(b)  Johnson-King 


(c)  RNG 


Figure  92.     Instantaneous  Particle  Trace,  Oscillatory  Motion,  M=0.3, 
<x=15°   (Downstroke),  a(t)=ll°±5°sin((0t),  k=0.10,  Re=4xl()6 
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EFFECT  OF  TURBULENCE  MODELING  ON  THE  COMPUTED 
TRAILING  EDGE  FLOWFIELD 

Moo  =  0.3,  a  (t)  =  11°  +  5°  sin  (cot),  k  =  0.1,  Re  =  4  x  106 
a  =  15^  DOWNSTROKE 


(a)  Baldwin-Lomax 


(b)  Johnson-King 


(c)  RNG 

Figure  93.     Velocity  Field,  Oscillatory  Motion,  M=0.3,  a=15 
(Downstroke),   a(t)=ll°±5°sin(a>t),  k=0.10,  Re=4xl()6 
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EFFECT  OF  TURBULENCE  MODELING  ON  THE  COMPUTED 
TRAILING  EDGE  FLOWFIELD 

M  =  0.3,  a  (t)  =  11°  +  5°  sin  (wt),  k  =  0.1,  Re  =  4  x  106 
a  =  14DOWNSTROKE 


(a)  Baldwin-Lomax 


(b)  Johnson-King 


(c)  RNG 


Figure  94.     Instantaneous  Particle  Trace,  Oscillatory  Motion,  M=0.3, 
<x=14°   (Downstroke),  a(t)=ll°±5°sin(o)t),  k=0.10,  Re=4xl06 
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EFFECT  OF  TURBULENCE  MODELING  ON  THE  COMPUTED 
TRAILING  EDGE  FLOWFIELD 

M  =  0.3,  a  (t)  =  11°  +  5°  sin  (cot),  k  =  0.1,  Re  =  4  x  106 

a  =  14°DOWNSTROKE 


(a)  Baldwin-Lomax 


(b)  Johnson-King 


(c)  RNG 

Figure  95.     Velocity  Field,  Oscillatory  Motion,  M=0.3,  a=14 
(Downstroke),   a(t)=ll°±5°sin(cot),  k=0.10,   Re=4xl()6 
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V.  CONCLUSIONS    AND    RECOMMENDATIONS 

In  the  preceding  investigation  a  computationally  efficient,  fully  factorized, 
two-dimensional  Navier-Stokes  flow  solver  was  utilized  to  predict  steady  and 
unsteady  flow  solutions  about  a  NACA  0012  airfoil.  Comparisons  between  a 
steady  viscous/inviscid  interaction  method  code  using  the  Cebeci-Smith 
turbulence  model  and  the  Navier-Stokes  code  using  the  Baldwin-Lomax 
turbulence  model  show  close  agreement  up  to  just  prior  to  the  static  stall  angle 
of  attack.  After  this  point,  the  viscous/inviscid  interaction  method  code  fails  to 
accurately  model  the  flow  separation. 

When  applied  to  the  unsteady  problem  of  airfoil  stall  flutter  in  compressible 
flow,  the  Navier-Stokes  code  shows  that  the  modeling  of  the  recirculatory  flow 
region  with  current  turbulence  models  fails  to  capture  the  essential  physics 
which  governs  the  onset  of  stall  flutter.  Comparisons  of  the  computed  results 
with  available  experimental  data  indicates  that  even  though  the  lift  response  is 
fairly  well  predicted,  the  computation  of  the  pitching  moment  hysteresis  loops  is 
very  sensitive  to  turbulence  modeling.  Results  computed  with  several  current 
models  are  in  good  agreement  whenever  the  steady  stall  angle  is  exceeded  only 
slightly.  However,  they  fail  to  capture  a  vortex  shedding  process  that  may 
contribute  to  the  onset  of  stall  flutter. 

Therefore,  further  detailed  studies  of  improved  numerical  schemes  and 
turbulence  models  as  well  as  viscous/inviscid  interaction  approaches  are 
required  to  improve  the  prediction  of  the  unsteady  flow  separation  and  vortex 
shedding  phenomenon.  Further  computational  work  with  the  full  Navier-Stokes 
equations  instead  of  applying  the  thin-layer  approximation  is  recommended.    In 
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addition,  local  grid  refinement  studies  that  focus  on  the  the  leading  and  trailing 
edge  upper  surface  may  lead  to  the  accurate  prediction  of  the  extremely 
complex  vortex  development  and  shedding  process. 
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APPENDIX   A   -  FLUX  JACOBIAN  MATRIX 

The  Jacobian  matrices  A=  oE/  aq  and  B  =  aF/  dq  are  given  by  A  or  B  = 
kol  +  kiA  +  k2B,  where  A  =  3E/3q  and  B  =  3F/3q  are  the  usual  Jacobian 
matrices  of  the  Cartesian  flux  vectors.   The  A  or  B  matrix  is 


r 


A 

dq 


k  k  k 

fk.0  IVj  A-2 

-uCk^+kjV)  -(y-2)k,u  -(y-ljkjv 

+  kj02  +  k0  +  k^*  k2v  +  k2u 


-v(kjU  +k2v) 
+  k2<t>2 


k,v 


-(y-l)k2u 


-(y-2)k2v  + 
k0  +kjU+k2v 


0         1 


(y-ijk, 


(y-i)k: 


(klU+k2v)*  [y(e/p)_02|kl-      [r(e/p)-02jk2-     y(klU  +  k2v) 

[[-y(e/p)+2^]    (y-i)(klU  +  k2v)u    (y-l)(klU  +  k2v)v    +  k0 


where  *2  =  0.5(y  -  l)(u2  +  v2),  ko=£t,  ki=(U,  k2=^z  for  A,  and  ko=Ct,  kl=Cx, 
k2<zforB. 
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APPENDIX  B   -  PITCHING   MOMENT  DERIVATION 

Figure  96  shows  the  nomenclature  for  the  integration  of  the  pressure  and 
shear  stress  distributions  over  a  two  dimensional  body  surface  as  generated  by 
the  Navier-Stokes  code  Ns2.f.  For  unsteady  motion  the  airfoil  is  rotated  with 
respect  to  a  fixed  inertial  frame  of  reference  (x,z:  Laboratory  frame  of 
Reference).  Note,  for  unsteady  motion,  as  the  airfoil  rotates  through  an  angle  of 
attack,  all  grid  (i=l...Imax>  k=l...Kmax)  are  rotated  about  the  desired  pivot  point 
with  respect  to  the  fixed  x,z  coordinate  system.  As  the  solution  is  advanced  in 
time  the  flow  quantities  for  this  new  grid  location  are  updated  for  each  point. 


s 


s 


Az=z(i+l)-z(i)f^^      *         H" 
Ax=x(i+1)  -  x(i) 


T.E.  lower,  i=3 1 


x  -  coordinate 
Figure  96.     Pitching  Moment  Derivation 
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PITCHING   MOMENT   DERIVATION   GENERATED   BY   MATHCAD 

1 .  Read  tabulated  data  generated  from  Ns2.f.  Data  computed  at  cc=0°  and  a  Re=4x  10**6. 


Data    :=  READPRN(up)     Data    :=  READPRN(low)       i  :=  1 
u  1 


50 


<0> 
x    :=  Data 
u  u 

<1> 
z    :=  Data 
u  u 

<2> 
cp    :=  Data 
u  u 

<3> 
t    :=  Data 
u  u 

2.  Compute  Ax  and  Az  for  upper  and  lower  surface. 


<0> 


x  :=  Data 

1  1 

z  :=  Data 

1  1 

cp  :=  Data 

1  1 

t  : =  Data 

1  1 


<1> 


<2> 


<3> 


Ax     :  =  x     -  x 

u        u       u 

i        i       i-1 


A  x  .  1   :=x.l   -x.l 

i        i       (i-D 

Az     :=z     -z  Az.l:=z.l-z.l 

u        u       u 

i         i       i-1  i         i        (i-1) 

3.  Compute  segment  length  As  for  upper  and  lower  surface. 


As 


u 


Ax 


u 


Az 


As.l 


2  2 

"A  x  .  1    ]        +     |"A  z  .  1 

>    .L  iJ  L  iJ    . 


4.  Compute  angle,  0,  relative  to  perpendicular  for  upper  and  lower  surface. 


Az 


u 


=  0 


+  atan 


=  0 


+  atan 


180 


Ax 


180 


Az 


Ax 


u 
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Computed  the  x  and  z  moment  arms  for  the  upper  and  lower  surface. 


X  +      X 

u  u 

i  i-1 


Z  +      Z 

U  U 

i  i-1 


XX 


zz 


u 


u 


X  +     X 

1  1 

i  i-1 


z  +    z 

1  1 

i  i-1 


XX 


zz 


.  Computepitching  moment  coefficient  about  quarter  chord  by  applying  equation  1.11  from  [Ref. 
,p.  17]. 


Cm 


25 


■Z 


-cp     "COS 
u 


cp    ' sin 
u 
i 

cp    *  cos 

1 


u 


u 


-cp    "sin 
1 


+  t    *  sin 
u 
i 

+  t    *  cos 
u 
i 

+  t    'sin 

1 
i 


+  T      'COS 
1 

i 


XX 


u 


1  ... 


-As 


u 


zz 


u 


XX 


•As 


zz 


Cm     =  0.005 
25 

.  The  computed  pitching  moment  coefficient  from  Ns2d.f  (Computed  on  a  Cray  Super-Computer) 
/as  Cm=0.0006,  as  compared  to  the  calculations  made  on  a  Macintosh  SE/30  where  Cm  =0.005. 
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Upper  Surface  Data: 

X 

z 

-Cp 

Tau 

-2.5e-l 

-1.6e-3 

-1.0319e+Q 

3e-3 

-2.497e-l 

3e-3   -9. 

841e-l        -7.4e- 

-3 

-2.479e-l 

7.6e-3 

-6.798e-l 

-1.5e-2 

-2.447e-l 

1.2e-2 

-4.091e-l 

-1.64e-2 

-2.404e-l 

1.65e-2 

-1.766e-l 

-1.54e-2 

-2.346e-l 

2.09e-2 

9.93e-2 

-1.19e-2 

-2.273e-l 

2.49e-2 

2.286e-l 

-8.4e-3 

-2.187e-l 

2.88e-2 

2.769e-l 

-6.3e-3 

-2.088e-l 

3.26e-2 

3.316e-l 

-4.9e-3 

-1.975e-l 

3.61e-2 

3.715e-l 

-3.8e-3 

-1.85e-l 

3.94e-2 

3.951e-l 

-3e-3 

-1.713e-l 

4.26e-2 

4.105e-l 

-2.4e-3 

-1.564e-l 

4.55e-2 

4.199e-l 

-1.9e-3 

-1.403e-l 

4.81e-2 

4.245e-l 

-1.5e-3 

-1.23e-l 

5.05e-2 

4.248e-l 

-1.2e-3 

-1.047e-l 

5.26e-2 

4.221e-l 

-le-3 

-8.54e-2 

5.44e-2 

4.168e-l 

-7e-4 

-6.51e-2 

5.6e-2 

4.094e-l 

-6e-4 

-4.38e-2 

5.73e-2 

3.996e-l 

-4e-4 

-2.17e-2 

5.82e-2 

3.904e-l 

-3e-4 

1.3e-3 

5.89e-2 

3.778e-l 

-2e-4 

2.5e-2 

5.93e-2 

3.628e-l 

-le-4 

4.95e-2 

5.94e-2 

3.494e-l 

Oe+0 

7.45e-2 

5.92e-2 

3.343e-l 

le-4 

1.002e-l 

5.88e-2 

3.184e-l 

le-4 

1.264e-l 

5.81e-2 

3.031e-l 

2e-4 

1.53e-l 

5.72e-2 

2.875e-l 

2e-4 

1.8e-l 

5.6e-2 

2.71e-l 

2e-4 

2.073e-l 

5.46e-2 

2.544e-l 

3e-4 

2.348e-l 

5.3e-2 

2.378e-l 

3e-4 

2.626e-l 

5.13e-2 

2.211e-l 

3e-4 

2.904e-l 

4.93e-2 

2.045e-l 

3e-4 

3.182e-l 

4.72e-2 

1.884e-l 

3e-4 

3.461e-l 

4.5e-2 

1.716e-l 

4e-4 

3.738e-l 

4.27e-2 

1.559e-l 

4e-4 

4.014e-l 

4.02e-2 

1.421e-l 

4e-4 

4.287e-l 

3.77e-2 

1.252e-l 

4e-4 

4.557e-l 

3.5e-2 

1.067e-l 

4e-4 

4.823e-l 

3.23e-2 

9.05e-2 

4e-4 

5.085e-l 

2.96e-2 

7.46e-2 

4e-4 

5.342e-l 

2.68e-2 

5.68e-2 

4e-4 

5.593e-l 

2.4e-2 

3.86e-2 

4e-4 

5.838e-l 

2.12e-2 

2.06e-2 

4e-4 

6.076e-l 

1.84e-2 

1.3e-3 

3e-4 

6.306e-l 

1.56e-2 

-2.05e-2 

3e-4 

6.528e-l 

1.29e-2 

-4.42e-2 

3e-4 

6.742e-l 

1.01e-2 

-7.02e-2 

3e-4 

6.946e-l 

7.4e-3 

-9.2e-2 

2e-4 

7.141e-l 

4.7e-3 

-1.441e-l 

le-4 

7.325e-l 

2.2e-3 

-1.088e-l 

Oe+0 

7.499e-l 

Oe+0   -2. 

34e-l            Oe+0 
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Lower  Surface  Data: 

x z -Cp Tau 


-2.5e-l 

-1.6e-3 

-1.0319e+0 

3e-3 

-2.486e-l 

-6.2e-3 

-7.832e-l 

4.9e- 

■3 

-2.458e-l 

-1.06e-2 

-4.815e-l 

2.3e- 

■3 

-2.419e-l 

-1.51e-2 

-2.554e-l 

-2e-4 

-2.366e-l 

-1.96e-2 

2.4e-2 

-2.2e- 

•3 

-2.298e-l 

-2.37e-2 

2.063e-l 

-2.8e- 

■3 

-2.215e-l 

-2.76e-2 

2.624e-l 

-2.7e- 

■3 

-2.12e-l 

-3.14e-2 

3.142e-l 

-2.5e- 

■3 

-2.012e-l 

-3.5e-2 

3.619e-l 

-2.2e- 

•3 

-1.891e-l 

-3.84e-2 

3.892e-l 

-1.9e- 

•3 

-1.758e-l 

-4.16e-2 

4.067e-l 

-1.7e- 

■3 

-1.613e-l 

-4.45e-2 

4.178e-l 

-1.4e- 

■3 

-1.456e-l 

-4.73e-2 

4.239e-l 

-1.2e- 

•3 

-1.288e-l 

-4.97e-2 

4.256e-l 

-le-3 

-1.109e-l 

-5.19e-2 

4.237e-l 

-8e-4 

-9.19e-2 

-5.39e-2 

4.191e-l 

-7e-4 

-7.2e-2 

-5.55e-2 

4.126e-l 

-5e-4 

-5.11e-2 

-5.69e-2 

4.035e-l 

-4e-4 

-2.94e-2 

-5.79e-2 

3.932e-l 

-3e-4 

-6.8e-3 

-5.87e-2 

3.84e-l 

-2e-4 

1.66e-2 

-5.92e-2 

3.681e-l 

-le-4 

4.06e-2 

-5.94e-2 

3.527e-l 

Oe+0 

6.53e-2 

-5.93e-2 

3.419e-l 

le-4 

9.06e-2 

-5.9e-2 

3.249e-l 

le-4 

1.164e-l 

-5.84e-2 

3.084e-l 

2e-4 

1.427e-l 

-5.76e-2 

2.941e-l 

2e-4 

1.694e-l 

-5.65e-2 

2.777e-l 

3e-4 

1.963e-l 

-5.52e-2 

2.612e-l 

3e-4 

2.236e-l 

-5.37e-2 

2.448e-l 

3e-4 

2.51e-l 

-5.2e-2 

2.283e-l 

4e-4 

2.785e-l 

-5.02e-2 

2.118e-l 

4e-4 

3.061e-l 

-4.82e-2 

1.957e-l 

4e-4 

3.336e-l 

-4.6e-2 

1.794e-l 

4e-4 

3.61e-l 

-4.38e-2 

1.631e-l 

4e-4 

3.883e-l 

-4.14e-2 

1.4  9e-l 

4e-4 

4.154e-l 

-3.89e-2 

1.342e-l 

4e-4 

4.421e-l 

-3.64e-2 

1.16e-l 

4e-4 

4.685e-l 

-3.37e-2 

9.88e-2 

5e-4 

4.945e-l 

-3.11e-2 

8.37e-2 

5e-4 

5.199e-l 

-2.84e-2 

6.73e-2 

5e-4 

5.448e-l 

-2.57e-2 

4.95e-2 

5e-4 

5.69e-l 

-2.29e-2 

3.18e-2 

4e-4 

5.926e-l 

-2.02e-2 

1.43e-2 

4e-4 

6.155e-l 

-1.75e-2 

-5.3e-3 

4e-4 

6.375e-l 

-1.48e-2 

-2.69e-2 

4e-4 

6.587e-l 

-1.21e-2 

-5.06e-2 

4e-4 

6.789e-l 

-9.5e-3 

-7.47e-2 

4e-4 

6.982e-l 

-6.9e-3 

-9.64e-2 

3e-4 

7.165e-l 

-4.4e-3 

-1.491e-l 

2e-4 

7.338e-l 

-2.1e-3 

-1.101e-l 

le-4 

7.499e-l 

Oe+0   -2. 

34e-l             Oe+0 
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APPENDIX  C   -  USING  THE  PROGRAMS 


A.   DIRECTIONS  FOR  THE  EXECUTION  OF  NS2.F 

1.     Steady    State   Case 

In  order  to  generate  a  Cp,  Q,  or  Cm  plot  from  Ns2.f  there  are  many  steps,  using 
six  different  programs  that  must  be  taken:  (1.)  grape  (grid  generation  program),  (2.)  rotgr 
(grid  rotation  program),  (3.)  Ns2.f  (Navier-Stokes  code),  (4.)vi  (editor),  (5.)  plot3d 
(flow  visualization  program),  and  (6.)xypIot  (a  simple  x-y  plot  program). 

A.)  Ns2.f  is  a  fortran  code  and  will  take  hours  to  run  even  on  the  Stardent  mini- 
super  computer.  To  run  this  code,  an  executable  file  (ns2),  an  input  file  ns2.in  (actual 
name  of  input  file  is  users  choice)  and  a  grid  file  FOR01  l.DAT  (Stardent)  or  fort.  1 1  (on 
other  computers)  must  be  in  the  directory. 

B.)  To  obtain  a  grid,  a  grid  generation  program  such  as  grape  must  first  be 
used.  Grape  output  places  the  leading  edge  of  the  grid  at  the  reference  axis  at  0.0°  AOA. 
Now,  the  grid  must  be  rotated  to  the  desired  AOA  by  using  the  program  rotgr.  The  input 
file  (output  from  grape)  for  rotgr  must  be  FOR001.DAT  and  the  output  is  FOR002.DAT. 
The  output  from  rotgr  (FOR002.DAT)  must  be  renamed  to  FOR01  l.DAT  for  use  in  ns2. 

C.)  Next  edit  the  input  file  to  reflect  the  desired  flow  conditions.  Any  editor  will 
do,  however  vi  is  perhaps  the  most  common  editor  and  it  is  found  on  most  unix  based 
machines.  Enter  Mach  No.  and  AOA,  check  smoothing  factors,  set  time  step,  Courant  No. 
and  initial  no.  of  iterations  (small  value-  for  first  run).  Check  restart  false  for  first  run  and 
make  sure  No.  point  of  your  grid  match  input  values. 

D.)  The  code  is  executed  as  follows:  ns2    <ns2.in>    ns2.out 

The  first  ten  iterations  should  take  1-10  minutes,  depending  on  the  type  of 
computer  used.  Several  output  files  will  be  generated: 

TABLE   3.   NS2.F    OUPUT 


OUTPUT 

COMMENTS 

ns2.out: 

convergence  data 

FOR004.DAT 

Qdata 

FOR007.DAT 

residuals 

FOR008.DAT 

Q,  Cd,  Cm  data 

FOR009.DAT 

Cp  data  (to  be  used  with  xyplot) 

FOR031.DAT 

Flow  data  (to  be  used  with  plot3d) 

E.)  Re-run  the  code  for  a  greater  number  of  iterations  to  achieve  a 
converged  solution.  When  the  change  in  the  max  residual  (FOR007.DAT)  has  dropped 
10~2,  the  solution  has  converged.  Ensure  restart  is  set  to  true.  The  CFL  can  also  be 
adjusted  to  obtain  a  converged  solution  more  rapidly. 
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a.    Sample  Input  to  Ns2.f 


MACH, 

ALFAO, 

ALFA1, 

ALFARE, 

REDFRE, 

REYNOLDS 

0.300, 

0.00, 

0.0, 

0.00, 

0.00, 

2.70 

ED2X, 

ED2Y, 

ED4X, 

ED4Y, 

ED 

0.00, 

0.00, 

0.030, 

0.030, 

2.0 

DT, 

COUR, 

NITER, 

NEWTIT 

0.0020, 

2100, 

50, 

2 

RSTRT, 

OSCIL, 

RAMP, 

NPER, 

TSHIFT 

true, 

false, 

false, 

10000, 

-0.5 

TIMEAC, 

IMPLBC, 

EXPLBC, 

CIRCOR 

1, 

false, 

true, 

false 

vise, 

BLTM, 

JKTM, 

RNGTM 

true, 

true, 

false, 

false 

I  TEL, 

ITEU 

31, 

131 

IA1,  IA2,  IA3,  IA4,  IA5,  IA6,  IA7,  IA8,  IA9,  IA10 
1250,1350,1450,1550,1570,1590,1610,1630,1630,  1805 
false 
UNSTST  (If  true  Time  =  0 . , Set  TRUE  only  when  unsteady  motion  starts  from  steady 
steady  state,  TRUE  for  steady-state  ok,  but  for  unsteady  restarts  must  be  FALSE 
for  proper  recording  of  unsteady  motion) 

Mach  :  Free  stream  Mach  number 

AlfaO  :  Angle  of  attack,  also  mean  angle  of  attack  for  unsteady 

Alfal  :  Amplitude  of  Oscillatory  motion 

Alfare  :  Reference  angle 

Redfre  :  Reduced  frequency  k  =  omega  *  c  /  2U 

Reynolds  :  Reynolds  number  Re  =  cU/n 

ED2x     :    X-direction  2nd  order  explicit  smoothing  (  e2x  =0.00  subsonic, 

0.25  <  e2x  <  0.50  transonic) 

ED2z     :    Z-direction  2nd  order  explicit  smoothing  (  e2z  =0.00  subsonic, 

0.10  <  e2z  <  0.20  transonic) 

ED4x     :    X-direction  4nd  order  explicit  smoothing  (  e4x  =0.03  subsonic, 

e4x  =0.05  transonic) 

ED4z     :    Z-direction  4nd  order  explicit  smoothing  (  e4z  =0.03  subsonic, 

e4z  =0.05  transonic) 

ED       :   Scaling  of  Implicit  smoothing 

ISPEC    :    Spectral  radious  parameter 

Dt  :  Time  step 

Cour  :  Courant  number  Cu  =  dt  *  L_max 

Niter  :  Number  of  Iteration  in  this  run 

Newtit  :  Newton  subiteration  within  each  timestep 

RSTRT  :  Restart 

OSCIL  :  Oscillatory  motion  A(t)  =  A0  +  Al  *  sin  (  k  *  M  *  t  ) 

RAMP  :  Ramp  motion 

NPER  :  Number  of  time  steps  in  one  period  of  oscillation,  dt=T/Nper 

TSHIFT  :  Time  shift  in  radians  to  start  oscilation  for  any  a(t) 

TIMEAC   :   Time  accurate  Tacc=l  and  for  Jacobian  Scaled  Dt,  Tacc=0 
IMPLBC    :    Implicit  wall  BC  Treatment 
EXPLBC   :   Explicit  wall  BC  Treatment 
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vise 

TURBL 

JKTM 

RNGTM 


Only  laminar  viscosity 
Baldwin-Lomax  eddy  viscosity 
Johnson-King  eddy  viscosity 
RNG  eddy  viscosity 


ITEL,    ITEU 
IA1  etc. 


Lower  and  Upper  trailing  edge  I  locations 
Write  out  (IAn/100)  angles  of  attack  in  units 
61  -  70 


Read  grid  from  unit  fort. 11  and  the  flow  from  unit  f ort . 31 
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b.    Sample  Output  from  Ns2.f 

Grid  dimensions  :  161x64 
Mach  =    0.30000 

Ao    =    0.00000    Al  =    0.00000    k  =    0.01272 
Re    =  0.2700E+07    Dt  =    0.00200    Cu=2100 . 00000 

Dimpl  =    2.00000 

D2x  =    0.00000   D2z  =    0.00000 

D4x  =    0.03000   D4z  =    0.03000 

Restart  =  T 
Oscil  =  F 
Ramp    =     F 

Itel  =    31   Iteu  =  131   He  =   81 

Timeac  =  T 
BC  impl  =  F 
BC  expl   =      T 

Cour   -  2100. 

L_max  =  232096.4825513 

Dt      =    0.00904796 

It    drmax  dumax  dvmax         demax      ir   kr 

3000   0.573103E-07   0.371716E-07   0.563395E-07   0.129865E-06    140    63 
Ni=    1   0.903247E-08   0.711341E-08   0.803250E-08   0.219128E-07 

Alfa(t)  =    0.000  CI  =   0.007196  Cd  =   0.003260  Cm  =  -0.000692 
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2.      Unsteady    Case 

The  execution  of  an  unsteady  case  is  essentially  the  same  as  for  the 
steady-state  case. 

A.  First,  a  well  converged  steady-state  solution  must  be  obtained. 
This  solution  then  becomes  the  starting  point  for  the  unsteady  solution. 

B.  The  restart  is  left  set  to  true,  however,  to  reset  the  time  counter 
to  time=0,  the  UNSTST  variable  must  be  set  to  true  for  several  iterations  and 
then  back  to  false  for  normal  time  counting. 

C.  As  the  unsteady  motion  begins  to  march  in  time  Flow  data  and 
Grid  position  data  may  be  stored  at  user  determined  angles  of  attack.    To  store 
flow  data  at  a=12.5°,  a=13.5°,  a=14.5°,  etc...  enter  the  data  shown  in  Table  4. 
into  the  input  file. 


TABLE 

4.   ANGLE   OF   ATTACK 

INPUT 

IA1 

IA2 

IA3 

IA4 

IA5 

IA6 

IA7 

IA8 

IA9 

IA10 

1250 

1350 

1450 

1550 

1570 

1590 

1610 

1630 

1630 

1805 

D.  The  Flow  data  and  the  Grid  data  at  each  a  will  be  stored  in  data 
files  FOR021.DAT  through  FOR030.DAT.  The  Grid  and  Flow  data  can  be 
separated  using  a  simple  fortran  code  that  reads  the  data  and  then  writes  the 
data  to  two  separate  files.  Splgf.f  inputs  the  desired  file  ID  and  outputs  a  Grid 
file  and  Flow  file. 
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a.    Sample  Input  to  Ns2.f 

MACH, 

ALFAO, 

ALFA1, 

ALFARE, 

REDFRE, 

REYNOLDS 

0.300, 

0.00, 

15.54, 

0.00, 

0.01272, 

2.70 

ED2X, 

ED2Y, 

ED4X, 

ED4Y, 

ED 

0.00, 

0.00, 

0.030, 

0.030, 

2.0 

DT, 

COUR, 

NITER, 

NEWTIT 

0.0020, 

2100, 

50, 

1 

RSTRT, 

OSCIL, 

RAMP, 

NPER, 

TSHIFT 

true, 

false, 

true, 

3000, 

-0.5 

TIMEAC, 

IMPLBC, 

EXPLBC, 

CIRCOR 

1, 

false, 

true, 

false 

vise, 

BLTM, 

JKTM, 

RNGTM 

true, 

true, 

false, 

false 

ITEL, 

ITEU 

31, 

131 

IA1, 

IA2,  ia: 

,  IA4, 

EA5,  IA6, 

IA7, 

IA8, 

IA9,  IA 

1250,1350,1450,1550,1570,1590,1610,1630,1630,  1805 
false 
UNSTST  (If  true  Time  =  0 . , Set  TRUE  only  when  unsteady  motion  starts  from  steady 
steady  state,  TRUE  for  steady-state  ok,  but  for  unsteady  restarts  must  be  FALSE 
for  proper  recording  of  unsteady  motion) 


Mach 

AlfaO 

Alfal 

Alfare 

Redfre 

Reynolds 

ED2x 

ED2z 

ED4x 

ED4z 

ED 
ISPEC 


Free  stream  Mach  number 

Angle  of  attack,  also  mean  angle  of  attack  for  unsteady 

Amplitude  of  Oscillatory  motion 

Reference  angle 

Reduced  frequency  k  -   omege  *  c  /  2U 

Reynolds  number  Re  =  cU/n 

X-direction  2nd  order  explicit  smoothing  (  e2x  =0.00  subsonic, 

0.25  <  e2x  <  0.50  transonic) 

Z-direction  2nd  order  explicit  smoothing  (  e2z  =0.00  subsonic, 

0.10  <  e2z  <  0.20  transonic) 

X-direction  4nd  order  explicit  smoothing  (  e4x  =0.03  subsonic, 

e4x  =0.05  transonic) 

Z-direction  4nd  order  explicit  smoothing  (  e4z  =0.03  subsonic, 

Scaling  of  Implicit  smoothing 

Spectral  radious  parameter 


Dt 

Cour 
Niter 
Newtit 

RSTRT 

OSCIL 

RAMP 

NPER 

TSHIFT 

TIMEAC 
IMPLBC 
EXPLBC 


Time  step 

Cour ant  number  Cu  =  dt  *  L_max 

Number  of  Iteration  in  this  run 

Newton  subiteration  within  each  timestep 

Restart 

Oscillatory  motion  A(t)  =  A0  +  Al  *  sin  (  k  *  M  *  t  ) 

Ramp  motion 

Number  of  time  steps  in  one  period  of  oscillation,  dt=T/Nper 

Time  shift  in  radians  to  start  oscilation  for  any  a(t) 

Time  accurate  Tacc=l  and  for  Jacobian  Scaled  Dt,  Tacc=0 
Implicit  wall  BC  Treatment 
Explicit  wall  BC  Treatment 
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vise 

TURBL 

JKTM 

RNGTM 


Only  laminar  viscosity 
Baldwin-Lomax  eddy  viscosity 
Johnson-King  eddy  viscosity 
RNG  eddy  viscosity 


ITEL,    ITEU 
IA1  etc. 


Lower  and  Upper  trailing  edge  I  locations 
Write  out  (IAn/100)  angles  of  attack  in  units 
61  -  70 


Read  grid  from  unit  fort. 11  and  the  flow  from  unit  fort. 31 
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b.    Sample   Output  from  Ns2.f 


Grid  dimensions  :  161x64 
Mach  =    0.30000 


Ao 
Re 

0.00000    Al 
0.2700E+07     Dt 

16.00000 
0.01829 

k  = 
Cu= 

0.01272 
15.00000 

Dimpl  = 
D2x  = 
D4x  = 

2.00000 
0.00000   D2z  = 
0.03000   D4z  = 

0.00000 
0.03000 

Restart 

Oscil 

Ramp 

T 
F 
T 

Itel  =    31  Iteu  =  131  He  =   81 

Timeac  =  T 
BC  impl  =  F 
BC  expl   =      T 

Cour  =  4289.735274387 
L_max  =  234539.9275225 
Dt     =    0.01829000 

It    drmax  dumax  dvmax         demax      ir   kr 

1950   0.849290E-06   0.107934E-05   0.190543E-05   0.174200E-05    119    63 
Ni=    1   0.113997E-06   0.302620E-06   0.448250E-06   0.196797E-06 

it  =  1950  Time  =  35.665500  omega  =  0.007632  phase  =   0.043322 
Alfa(t)  =   15.596  CI  -   1.425268  Cd  -    0.137214  Cm  =    0.003384 

2000   0.778403E-06   0.107038E-05   0.191827E-05   0.158330E-05     118    63 
Ni=    1   0.113176E-06   0.297453E-06   0.450476E-06   0.197668E-06 

it  =  2000  Time  =  36.580000  omega  =   0.007632  phase  =   0.044433 
Alfa(t)  =   15.996  CI  =    1.453050  Cd  =    0.146989  Cm  =    0.001525 
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B.    DIRECTIONS  FOR  THE  EXECUTION  OF  U2DIIF.F 
1.      Unsteady    Case 

A.    A  complete  description  of  the  input  and  output  parameters  can  be 
found  in  Teng. 

TABLE   5.   U2DIIF.F    INPUTS/OUTPUTS 


FILE 

COMMENTS 

FOR001.DAT 

Input 

UMif.out 

Output 

B.)  For  a  unix  based  machine  the  code  is  executed  as  follows: 

U2diif  >  U2diif.out 

Note:  The  original  Fortran  code  can  be  modified  to  store  desired  output  in 
separate  data  files  by  inserting  a  WRITE  (9,*)  statement  inside  the  required  do- 
loop.  (Example:  Cp  vs  x/c  data  could  be  stored  in  a  data  file  FOR009.DAT). 
This  requires  some  knowledge  of  Fortran  programming. 
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a.    Sample  Input  to  U2diif.f 

4 

AIRFOIL  TYPE  :  NACA  0012  AIRFOIL 

NLOWER  =  50  ,  NUPPER  =  50 
******************************************************* 

0,  50,  50 

12 

0.00,  15.54,  10.67,  0.00,  0.25,  0.00,  0.00 

0.00,  0.00,  0.00 

11.00,  0.010,  0.0001,  0.00 

108,  294,  487,  584,  802,  891,  1006,  1007,  1283,  1554 

ITITLE 

INPUT  TITLE  (ITITLE  =  NO.  LINE) 

IFLAG   NLOWER   NUPPER 

AIRFOIL  TYPE 

ALPI   DALP   TCON   FREQ   PIVOT   UGUST   VGUST 

DELHX   DELHY   PHASE 

TF   DTS   TOL   TADK 

ial  ia2  ia3  ia4  ia5  ia6  ia7  ia8  ia9  ialO 
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b.     Sample  Output  to  U2diif.f 


DATA  FROM  f ort . 1 : (IRIS)  OR  FOR001.DAT: (STARDENT) 

****************************************************** 

AIRFOIL  TYPE  :  NACA  0012  AIRFOIL 
NLOWER  =  50  ,  NUPPER  =  50 
****************************************************************** 


IFLAG  (0:NACA,  1: INPUT)  =  0 
NO.  PANELS  UPPER  SURFACE  =  50 
NO.  PANELS  LOWER  SURFACE      =    50 


INITIAL  ANGLE  OF  ATTACK  =  0.0000 

FINAL    ANGLE  OF  ATTACK  =  15.5400 

RISE  TIME  =  10.6700 

REDUCED  FREQ.  FOR  OSCIL  =  0.0000 

PIVOT  POINT  =  0.2500 

STREAMWISE  GUST  VELOCITY  =  0.0000 

PERPENDICULAR  GUST  VELOCITY   =  0.0000 


X  AMPLITUDE  OF  TRANS  OSCILL.  =  0.0000 
Y  AMPLITUDE  OF  TRANS  OSCILL.  =  0.0000 
PHASE        OF  TRANS  OSCILL.  =     0.0000 


FINAL  TIME                    =  11.0000 

INITIAL  TIME  STEP  FOR  S.S.    =  0.0100 

TOLERANCE  FOR  CONVERGENCE     =  0.0001 

FACTOR  BY  WHICH  DTS  ADJUSTED  =  0.0000 


COORDINATES  OF  AIRFOIL  NODES 

X/C        Y/C 
1.000000   0.000000 
0.999013  -0.000141 
0.996057  -0.000562 


0.984292  0.002222 

0.991144  0.001258 

0.996057  0.000562 

0.999013  0.000141 

1.000000  0.000000 

AIRFOIL  PERIMETER  LENGTH  =    2.039290 

STEADY  FLOW  SOLUTION  AT  ALPHA  =    0.000000 

J     X(J)       Q(J)      GAMMA      CP(J)       V(J) 
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1 

2 
3 
4 
5 
6 
7 
8 
9 
10 


0.999507 
0.997535 
0.993600 
0.987718 
0.979910 
0.970208 
0.958651 
0.945283 
0.930159 
0.913336 


-0 


0.112073 
0.120842 
0.126091 
0.129448 
0.131624 
0.132948 
0.133600 
0.133699 
0.133328 
132558 


0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 


0.433845 
0.341896 
0.279513 
0.232223 
0.193468 
0.160181 
0.130705 
0.104047 
0.079564 


0.056812   - 


0.752433 
0.811236 
0.848815 
0.876229 
0.898071 
0.916416 
0.932360 
0.946548 
0.959393 
0.971179 


90 

0, 

.894883 

-0, 

.131445 

0, 

.000000 

0, 

,035461 

0. 

,982110 

91 

0, 

,913336 

-0, 

.132561 

0, 

.000000 

0, 

,056813 

0, 

,971178 

92 

0, 

.930159 

-0, 

.133333 

0, 

.000000 

0, 

,079566 

0, 

,959393 

93 

0, 

.945283 

-0, 

.133704 

0, 

.000000 

0, 

,104048 

0, 

,946547 

94 

0. 

.958651 

-0, 

.133607 

0, 

.000000 

0, 

,130706 

0, 

,932359 

95 

0, 

,970208 

-0. 

.132956 

0, 

,000000 

0, 

,160183 

0. 

,916415 

96 

0, 

.979910 

-0, 

,131634 

0, 

,000000 

0, 

,193471 

0, 

,898070 

97 

0. 

,987718 

-0, 

,129461 

0, 

,000000 

0, 

,232226 

0. 

,876227 

98 

0. 

,993600 

-0. 

,126109 

0, 

,000000 

0. 

,279518 

0. 

,848812 

99 

0. 

,997535 

-0, 

,120869 

0, 

,000000 

0, 

,341902 

0. 

,811232 

100 

0. 

999507 

-0. 

,112063 

0. 

,000000 

0. 

,433845 

0, 

,752433 

COMPARISON 


0  F 


GAMMAS 


GAMMA  FROM  KUTTA  CONDITION: 
GAMMA  FROM  CONTOUR  INTEGR  (TRAIL  EDGE) : 
GAMMA  FROM  CONTOUR  INTEGR  (MIDPOINTS) : 
GAMMA  FROM  BOX  INTEGR  (OFF  THE  CONTOUR) 
GAMMA  FROM  PRECISE  CONTOUR  INTEG  (6  PT) 


-0.00000024 
-0.00000025 
-0.00000025 
-0.00000024 
-0.00000024 


CD  =   0.000324 


CL  =  -0.000001 


CM  =   0.000001 


•  ••♦•••••••it*************************** 

***   BEGIN  UNSTEADY  FLOW  SOLUTION   **** 
••••••a******************************** 


TIME  STEP  TK  = 


0.010000 


TK  -  TKM1  = 


0.010000 


ALPHA (T)  = 
U(T) 


NITR 
0 
1 
2 
3 


0.000041 
0.000000 


VXW 
1.000000 
0.834730 
0.827773 
0.827452 


VYW 
0.000000 
0.000244 
0.000255 
0.000256 


OMEGA (T) 
V(T) 

WAKE 
0.010000 
0.008347 
0.008278 
0.008275 


=  -0.000143 
=  -0.000036 

THETA      GAMMA 
0.000000  -0.240067E-06 
0.000292   0.436698E-05 


0.000309 
0.000309 


0.412307E-05 
0.411209E-05 


CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  = 
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J 

X(J) 

Q(J) 

GAMMA 

CP(J) 

V(J) 

PREV  PHI (J) 

CURR  PHI  (J 

1 

0 

999507 

-0 

107774 

0 

.000004 

0.432728 

-0 

753023 

0 

032660 

0 

032661 

2 

0 

997535 

-0 

117297 

0 

.000004 

0.341094 

-0 

811709 

0 

033045 

0 

033045 

3 

0 

993600 

-0 

123077 

0 

000004 

0.279152 

-0 

849194 

0 

033663 

0 

033661 

4 

0 

987718 

-0 

126843 

0 

000004 

0.232351 

-0 

876535 

0 

034408 

0 

034404 

5 

0 

979910 

-0 

129345 

0 

000004 

0.194095 

-0 

898321 

0 

035211 

0 

035205 

6 

0 

970208 

-0 

130932 

0 

000004 

0.161293 

-0 

916622 

0 

036018 

0 

036011 

7 

0 

958651 

-0 

131801 

0 

000004 

0.132278 

-0 

932531 

0 

036785 

0 

036776 

8 

0 

945283 

-0 

132080 

0 

000004 

0.106051 

-0 

946691 

0 

037475 

0 

037463 

9 

0 

930159 

-0 

131864 

0 

000004 

0.081967 

-0 

959514 

0 

038055 

0 

038041 

0 

0 

913336 

-0 

131225 

0 

000004 

0.059576 

-0 

971280 

0 

038497 

0 

038482 

90 

0 

894883 

-0 

132662 

0 

000004 

0 

032373 

0 

982025 

0 

038778 

0 

038794 

91 

0 

913336 

-0 

.133893 

0 

000004 

0 

054048 

0 

971077 

0 

038497 

0 

038511 

92 

0 

930159 

-0 

134798 

0 

000004 

0 

077163 

0 

959272 

0 

038054 

0 

038067 

93 

0 

945283 

-0 

135322 

0 

000004 

0 

102044 

0 

946404 

0 

037474 

0 

037486 

94 

0 

958651 

-0 

135406 

0 

000004 

0 

129134 

0 

932188 

0 

036785 

0 

036794 

95 

0 

970208 

-0 

134972 

0 

000004 

0 

159071 

0 

916209 

0 

036018 

0 

036025 

96 

0 

979910 

-0 

133914 

0 

000004 

0 

192844 

0 

897820 

0 

035210 

0 

035216 

97 

0 

987718 

-0 

132066 

0 

000004 

0 

232097 

0 

875921 

0 

034407 

0 

034411 

98 

0 

993600 

-0 

129123 

0 

000004 

0 

279879 

0 

848433 

0 

033662 

0 

033664 

99 

0 

997535 

-0 

124414 

0 

000004 

0 

342703 

0 

810760 

0 

033045 

0 

033044 

100 

0. 

999507 

-0 

116362 

0 

000004 

0 

434960 

0 

751843 

0 

032659 

0 

032658 

COMPARISON    OF    GAMMAS 

GAMMA  FROM  KUTTA  CONDITION:  0.00000411 

GAMMA  FROM  CONTOUR  INTEGR  (TRAIL  EDGE):  -0.00000165 

GAMMA  FROM  CONTOUR  INTEGR  (MIDPOINTS):  -0.00000137 

GAMMA  FROM  BOX  INTEGR  (OFF  THE  CONTOUR):  0.00000462 

GAMMA  FROM  PRECISE  CONTOUR  INTEG  (6  PT) :  0.00000385 


CD  =   0.000324 


CL 


0.005136 


CM 


-0.003050 


TRAILING  VORTICES  DATA 

M    X(M)       Y(M)       CIRC 
1   1.004137   0.000001  -0.000009 


TIME  STEP  TK 


0.020000 


TK  -  TKM1  =    0.010000 


ALPHA(T)  =    0.000164 
U(T)      =    0.000000 


OMEGA (T)  =   -0.000285 
V(T)      =   -0.000071 


NITR     VXW        VYW       WAKE       THETA      GAMMA 

0   0.827452   0.000256   0.008275   0.000309   0.411158E-05 
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CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  =    0 

J     X(J)       Q(J)      GAMMA      CP(J)       V(J)  PREV  PHI  ( J)  CURR  PHI  (J) 

1  0.999507  -0.104926   0.000011   0.432928  -0.753413   0.032661  0.032658 

2  0.997535  -0.114799   0.000011   0.341275  -0.812045   0.033045  0.033042 


20  0.669285  -0.109493  0.000011  -0.121525  -1.061102  0.031441  0.031419 

21  0.639427  -0.105975  0.000011  -0.138604  -1.069082  0.029361  0.029339 

22  0.609018  -0.102029  0.000011  -0.155875  -1.077078  0.027013  0.026992 

23  0.578179  -0.097571  0.000011  -0.173377  -1.085108  0.024398  0.024378 
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C.    DIRECTIONS  FOR  THE  EXECUTION  OF  INCOMPBL.F 
1.     Steady    State   Case 

A.    Execution  of  Incompbl.f  is  straight  forward.    First  edit  the  input  files, 
FOR001.DAT  and  incompbl.dat,  and  set  the  desired  flow  conditions. 

TABLE  6.      INCOMPBL.F    INPUTS 


INPUT 

COMMENTS 

incompbl.dat 

IWAKE-Flag  Wake  Calulation 

NXT 

NW-No.  of  point  in  wake 

ITREND-No.  of  Iterations 

ITR(l)-Transition  Flag  Upper  Surface 

ITR(2)-Transition  Flag  Upper  Surface 

ISWPMX 

RL-Reynolds  Number 

XCTR(l)-Transition    Location    Upper 

Surface 

IP-Print  flag 

FOR001.DAT 

Input    AOA,    pivot    point    and    Airfoil 
Coordinates,  and  No.  of  Panels 

B.)  For  a  unix  based  machine  the  code  is  executed  as  follows: 
incompbl    <incompbI.dat>   incompbl.out. 
TABLE  7.      INCOMPBL.F   OUPUT 


OUTPUT 

COMMENTS 

incompbl.out 

Flow  data 

Note:  The  original  Fortran  code  can  be  modified  to  store  desired  output  in 
separate  data  files  by  inserting  a  WRITE  (24,*)  statement  inside  the  required  do-loop. 
(Example:  Cp  vs  x/c  data  could  be  stored  in  a  data  file  FOR024.DAT).  This  requires  some 
knowledge  of  Fortran  Programming. 
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a.    Sample  Input  to  Incompbl.f 
File  Named  fort.l  ofFOR001.DAT 

3 

NACA  0012  AIRFOIL 


ALP  I 

,000000 

NLOWER 

50 


X/C 


1.00000 
0.88000 
0.76000 
0.64000 
0.52000 
0.40000 
0.28000 
0.16000 
0.04000 
0.08000 
0.20000 
0.32000 
0.44000 
0.56000 
0.68000 
0.80000 
0.92000 
Y/C 

-0.00126 

-0.01694 

-0.03056 

-0.04222 

-0.05165 

-0.05803 

-0.05997 

-0.05444 

-0.03231 

0.04309 

0.05737 

0.05995 

0.05634 

0.04878 

0.03856 

0.02623 

0.01196 


PIVOT 

0.250000 

NUPPER 

50 

0.98000 
0.86000 
0.74000 
0.62000 
0.50000 
0.38000 
0.26000 
0.14000 
0.02000 
0.10000 
0.22000 
0.34000 
0.46000 
0.58000 
0.70000 
0.82000 
0.94000 

-0.00403 

-0.01935 

-0.03264 

-0.04396 

-0.05294 

-0.05868 

-0.05966 

-0.05236 

-0.02382 

0.04683 

0.05838 

0.05966 

0.05530 

0.04723 

0.03664 

0.02399 

0.00938 


0.96000 
0.84000 
0.72000 
0.60000 
0.48000 
0.36000 
0.24000 
0.12000 
0.00000 
0.12000 
0.24000 
0.36000 
0.48000 
0.60000 
0.72000 
0.84000 
0.96000 

-0.00674 

-0.02170 

-0.03467 

-0.04563 

-0.05415 

-0.05923 

-0.05911 

-0.04990 

0.00000 

0.04990 

0.05911 

0.05923 

0.05415 

0.04563 

0.03467 

0.02170 

0.00674 


File  named  incompbl.dat 


IWAKE 

0 

ITR(l) 

0 

IP 

1 


NXT 

161 

ITR(2) 

0 


NW 

37 

ISWPMX 

1 


0.94000 
0.82000 
0.70000 
0.58000 
0.46000 
0.34000 
0.22000 
0.10000 
0.02000 
0.14000 
0.26000 
0.38000 
0.50000 
0.62000 
0.74000 
0.86000 
0.98000 

-0.00938 

-0.02399 

-0.03664 

-0.04723 

-0.05530 

-0.05966 

-0.05838 

-0.04683 

0.02382 

0.05236 

0.05966 

0.05868 

0.05294 

0.04396 

0.03264 

0.01935 

0.00403 

I TREND 

40 

RL 

3000000.0 


0.92000 
0.80000 
0.68000 
0.56000 
0.44000 
0.32000 
0.20000 
0.08000 
0.04000 
0.16000 
0.28000 
0.40000 
0.52000 
0.64000 
0.76000 
0.88000 
1.00000 

-0.01196 

-0.02623 

-0.03856 

-0.04878 

-0.05634 

-0.05995 

-0.05737 

-0.04309 

0.03231 

0.05444 

0.05997 

0.05803 

0.05165 

0.04222 

0.03056 

0.01694 

0.00126 


XCTR(l) 
0.10000 


0.90000 
0.78000 
0.66000 
0.54000 
0.42000 
0.30000 
0.18000 
0.06000 
0.06000 
0.18000 
0.30000 
0.42000 
0.54000 
0.66000 
0.78000 
0.90000 


-0.01448 

-0.02842 

-0.04042 

-0.05026 

-0.05726 

-0.06006 

-0.05607 

-0.03842 

0.03842 

0.05607 

0.06006 

0.05726 

0.05026 

0.04042 

0.02842 

0.01448 
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b.    Sample   Output  from  Incompbl.f 

:     NACA  0012  AIRFOIL 
INPUT  DATA  FOR  INVISCID-FLOW  CALCULATIONS 


ALPI= 
NLOWER= 

COORDINATES 

x/c 

1.000000 
0.880000 
0.760000 
0.640000 
0.520000 
0.400000 
0.280000 
0.160000 
0.040000 
0.080000 
0.200000 
0.320000 
0.440000 
0.560000 
0.680000 
0.800000 
0.920000 

Y/C 

-0.001260 

-0.016940 

-0.030560 

-0.042220 

-0.051650 

-0.058030 

-0.059970 

-0.054440 

-0.032310 

0.043090 

0.057370 

0.059950 

0.056340 

0.048780 

0.038560 

0.026230 

0.011960 


4.0000   PIVOT= 
50   NUPPER- 


OF  THE  BODY 


0.980000 
0.860000 
0.740000 
0.620000 
0.500000 
0.380000 
0.260000 
0.140000 
0.020000 
0.100000 
0.220000 
0.340000 
0.460000 
0.580000 
0.700000 
0.820000 
0.940000 


0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 
0. 


960000 
840000 
720000 
600000 
480000 
360000 
240000 
120000 
000000 
120000 
240000 
360000 
480000 
600000 
720000 
840000 
960000 


0.2500 
50 


0.940000 
0.820000 
0.700000 
0.580000 
0.460000 
0.340000 
0.220000 
0.100000 
0.020000 
0.140000 
0.260000 
0.380000 
0.500000 
0.620000 
0.740000 
0.860000 
0.980000 


0.920000 
0.800000 
0.680000 
0.560000 
0.440000 
0.320000 
0.200000 
0.080000 
0.040000 
0.160000 
0.280000 
0.400000 
0.520000 
0.640000 
0.760000 
0.880000 
1.000000 


-0.004030 

■0.019350 

■0.032640 

■0.043960 

■0.052940 

-0.058680 

-0.059660 

-0.052360 

-0.023820 

0.046830 

0.058380 

0.059660 

0.055300 

0.047230 

0.036640 

0.023990 

0.009380 


-0.006740  -0. 

-0.021700  -0. 

-0.034670  -0. 

-0.045630  -0. 

-0.054150  -0. 

-0.059230  -0. 

-0.059110  -0. 

-0.049900  -0. 

0.000000   0. 

0.049900   0. 

0.059110   0. 

0.059230   0. 

0.054150   0. 

0.045630   0. 

0.034670   0. 

0.021700   0. 

0.006740   0. 


009380  -0.011960 

023990  -0.026230 

036640  -0.038560 

047230  -0.048780 

055300  -0.056340 

059660  -0.059950 

058380  -0.057370 

046830  -0.043090 

023820  0.032310 


052360 
059660 
058680 
052940 
043960 
032640 
019350 
004030 


0.054440 
0.059970 
0.058030 
0.051650 
0.042220 
0.030560 
0.016940 
0.001260 


0.900000 
0.780000 
0.660000 
0.540000 
0.420000 
0.300000 
0.180000 
0.060000 
0.060000 
0.180000 
0.300000 
0.420000 
0.540000 
0.660000 
0.780000 
0.900000 


-0.014480 

-0.028420 

-0.040420 

-0.050260 

-0.057260 

-0.060060 

-0.056070 

-0.038420 

0.038420 

0.056070 

0.060060 

0.057260 

0.050260 

0.040420 

0.028420 

0.014480 


INVISCID  FLOW  RESULTS 


PANEL 
1 
2 


XP 

0.99000E+00 
0.97000E+00 


YP 

-0.14950E-02 
-0.44813E-02 


CP 

0.25519E+00 
0.17805E+00 


174 


3 
4 
5 
6 
7 
8 
9 
10 


0.9500 
0.9300 
0.9100 
0.8900 
0.8700 
0.8500 
0.8300 
0.8100 


0E+00 
0E+00 
0E+00 
0E+00 
0E+00 
0E+00 
0E+00 
0E+00 


-0 


0.74415E-02 
0.10321E-01 
13061E-01 
0.15646E-01 
0.18113E-01 
0.20502E-01 
0.22829E-01 
0.25104E-01 


0.13179E+00 
0.97803E-01 
0.74658E-01 
0.60117E-01 
0.49465E-01 
0.39706E-01 
0.30916E-01 
0.21877E-01 


90 

0 

. 79000E+00 

0 

.  27319E-01 

-0, 

.12424E+00 

91 

0 

.81000E+00 

0 

.25096E-01 

-0, 

.10630E+00 

92 

0 

.83000E+00 

0 

.22823E-01 

-0. 

.  88121E-01 

93 

0 

. 85000E+00 

0 

.20492E-01 

-0, 

.  69791E-01 

94 

0 

.87000E+00 

0, 

.18084E-01 

-0, 

.49935E-01 

95 

0, 

.89000E+00 

0 

.15564E-01 

-0, 

.24583E-01 

96 

0, 

. 91000E+00 

0, 

.12899E-01 

0, 

,  93575E-02 

97 

0, 

, 93000E+00 

0, 

.10105E-01 

0, 

,  51153E-01 

98 

0, 

. 95000E+00 

0, 

.72366E-02 

0. 

,  99432E-01 

99 

0, 

, 97000E+00 

0, 

.43440E-02 

0. 

.15874E+00 

100 

0, 

. 99000E+00 

0, 

.14480E-02 

0. 

.25519E+00 

INVISCID  WAKE 

RESULTS 

PANEL 

XP 

YP 

CP 

101 

0, 

.10033E+01 

0. 

.26712E-04 

0. 

,  34783E+00 

102 

0. 

.10110E+01 

0, 

,  93657E-04 

0, 

.26481E+00 

103 

0, 

, 10209E+01 

0, 

.19595E-03 

0, 

.21845E+00 

104 

0, 

,10340E+01 

0, 

.35320E-03 

0. 

,18376E+00 

105 

0, 

.10510E+01 

0, 

,  59334E-03 

0, 

,15510E+00 

106 

0. 

.  10732E+01 

0. 

,  95766E-03 

0. 

,13035E+00 

107 

0. 

,11023E+01 

0, 

.15069E-02 

0, 

.10857E+00 

108 

0. 

.11402E+01 

0. 

.23302E-02 

0. 

,  89285E-01 

109 

0. 

.11897E+01 

0. 

.35565E-02 

0. 

.  72269E-01 

110 

0. 

12543E+01 

0. 

,  53713E-02 

0, 

.  57416E-01 

111 

0. 

13387E+01 

0, 

.80381E-02 

0. 

,44666E-01 

112 

0. 

.14489E+01 

0. 

,  11926E-01 

0. 

,33948E-01 

113 

0. 

15928E+01 

0, 

,  17545E-01 

0. 

.25162E-01 

114 

0. 

17806E+01 

0. 

.25591E-01 

0. 

.18168E-01 

115 

0. 

.20259E+01 

0. 

,  36999E-01 

0, 

,12767E-01 

116 

0. 

23462E+01 

0. 

,  53013E-01 

0. 

,87308E-02 

117 

0. 

27643E+01 

0. 

,  75273E-01 

0. 

,58084E-02 

118 

0. 

,  33103E+01 

0. 

.10592E+00 

0, 

37624E-02 

119 

0. 

38079E+01 

0. 

,13479E+00 

0. 

.26992E-02 

CL  =  0.47937E+00 


INPUT  DATA  FOR  BOUNDARY -LAYER  CALCULATIONS 


IWAKE 

NXT 

NW     ITREND 

0 

161 

37         40 

ITR(l) 

ITR(2) 

ISWPMX  10**-6*RL 

XCTR(l) 

0 

0 

1       3.00 

0.10 

IP 

1 

175 


******••*•****  CYCLE   40  ************** 

SUMMARY  OF  THE  DRAG,  LIFT  AND  PITCHING  MOMENT 
COEFFICIENTS  WITH  THE  CYCLE 

CD,CL  AND  CM  ARE  EVALUATED  FROM  THE  INTEGRATION  OF  CP 


CYCLE 

CD 

CL 

CM 

1 

0, 

,001978 

0, 

,001684 

0, 

,427162 

2 

0, 

.003184 

0. 

.001730 

0. 

,427756 

3 

0. 

.002734 

0. 

.001684 

0, 

,427162 

4 

0. 

,002422 

0. 

,001730 

0, 

,427770 

5 

0, 

.002227 

0. 

,001684 

0, 

,427157 

35 

0. 

,001684 

0, 

,427157 

0, 

.004733 

36 

0, 

,001731 

0, 

,427790 

0, 

.004494 

37 

0, 

,001684 

0, 

,427165 

0, 

.004731 

38 

0. 

,001731 

0, 

,427807 

0, 

,004491 

39 

0, 

,001684 

0, 

,427167 

0, 

,004731 

40 

0, 

,001731 

0. 

,427822 

0, 

,004488 

BOUNDARY  LAYER  PROPERTIES  FOR  THE  LAST  CYCLE 


UPPER  SURFACE  

XCTR=  0.134E+00 

NX        XC        XS        CF       DLS        UE        CP  IT 


76 

0, 

,00294 

0, 

,003411 

0, 

,06394 

0, 

,00005 

0, 

.15840 

0, 

,97491 

3 

77 

0, 

,00140 

0, 

,006562 

0, 

,03670 

0, 

,00005 

0, 

.33362 

0, 

,88870 

2 

78 

0, 

,00056 

0, 

,009021 

0, 

,02384 

0, 

,00005 

0. 

,48286 

0. 

,76685 

2 

79 

0, 

,00018 

0, 

,010828 

0. 

,01910 

0, 

,00005 

0. 

,59624 

0. 

,64450 

2 

158  0.98712  1.018435  0.00046  0.00711  0.87987  0.22583  5 

159  0.99209  1.023461  0.00035  0.00749  0.87320  0.23751  4 

160  0.99639  1.027802  0.00025  0.00785  0.86751  0.24742  5 

161  1.00000  1.031453  0.00019  0.00818  0.86283  0.25553  5 


LOWER  SURFACE  

XCTR=  0.677E+00 

NX        XC        XS        CF       DLS        UE        CP  IT 
87   0.00540  0.000436   0.47509   0.00006   0.02025   0.99959   2 
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88  0.00888  0.004971   0.15459   0.00006   0.18949   0.96409   2 

89  0.01333  0.010187   0.02621   0.00007   0.36027   0.87021   2 

90  0.01841  0.016077   0.01498   0.00007   0.53294   0.71597   3 


158  0.98713  0.992288  0.00204  0.00181  0.89064  0.20676  2 

159  0.99210  0.997313  0.00186  0.00191  0.88051  0.22471  2 

160  0.99639  1.001654  0.00171  0.00202  0.87121  0.24100  2 

161  1.00000  1.005305  0.00159  0.00212  0.86310  0.25505  2 
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